From ec972a7cdf7b4a529ad77a97cf178a7e2452619c Mon Sep 17 00:00:00 2001 From: tforest Date: Mon, 23 May 2022 20:25:12 +0200 Subject: [PATCH] update sfs plot and vcf parsing --- __init__.py | 3 +- bam_utils.py | 17 + customgraphics.py | 8 +- dependences/__init__.py | 1 + .../__pycache__/__init__.cpython-310.pyc | Bin 0 -> 222 bytes .../__pycache__/__init__.cpython-38.pyc | Bin 0 -> 223 bytes .../__pycache__/__init__.cpython-39.pyc | Bin 0 -> 192 bytes dependences/__pycache__/pybam.cpython-310.pyc | Bin 0 -> 38989 bytes dependences/__pycache__/pybam.cpython-38.pyc | Bin 0 -> 39196 bytes dependences/__pycache__/pybam.cpython-39.pyc | Bin 0 -> 39119 bytes dependences/pybam.py | 755 ++++++++++++++++++ dependences/pybam.py~ | 743 +++++++++++++++++ dependences/pybam_old.py | 743 +++++++++++++++++ sfs_tools.py | 73 +- stats_sfs.py | 1 - vcf_utils.py | 19 + 16 files changed, 2346 insertions(+), 17 deletions(-) create mode 100644 bam_utils.py create mode 100644 dependences/__init__.py create mode 100644 dependences/__pycache__/__init__.cpython-310.pyc create mode 100644 dependences/__pycache__/__init__.cpython-38.pyc create mode 100644 dependences/__pycache__/__init__.cpython-39.pyc create mode 100644 dependences/__pycache__/pybam.cpython-310.pyc create mode 100644 dependences/__pycache__/pybam.cpython-38.pyc create mode 100644 dependences/__pycache__/pybam.cpython-39.pyc create mode 100644 dependences/pybam.py create mode 100644 dependences/pybam.py~ create mode 100644 dependences/pybam_old.py diff --git a/__init__.py b/__init__.py index 65f0ffb..0f057a6 100644 --- a/__init__.py +++ b/__init__.py @@ -1 +1,2 @@ -from frst import sfs_tools, customgraphics, vcf_utils, sfs_tools, stats_sfs +from frst import sfs_tools, customgraphics, vcf_utils, sfs_tools, stats_sfs, dependences + diff --git a/bam_utils.py b/bam_utils.py new file mode 100644 index 0000000..7c2444c --- /dev/null +++ b/bam_utils.py @@ -0,0 +1,17 @@ +from frst.dependences import pybam + +def parse_bam(bam_file): + cov = {} + cov_val = 0 + for alignment in pybam.read(bam_file): + #pos 1 based (see Pybam doc) + start_pos = alignment.sam_pos1 + end_pos = alignment.sam_pos1 + alignment.sam_block_size + if start_pos not in cov: + cov[start_pos] = 0 + if end_pos not in cov: + cov[end_pos] = 0 + # add a weight on the start and lower on the end + cov[start_pos] = cov_val +1 + cov[end_pos] = cov_val -1 + return cov diff --git a/customgraphics.py b/customgraphics.py index 08bcd3f..60c833b 100644 --- a/customgraphics.py +++ b/customgraphics.py @@ -200,8 +200,12 @@ def scatter(x, y, ylab=None, xlab=None, title=None): plt.title(title) plt.show() -def barplot(x, y, ylab=None, xlab=None, title=None): - plt.bar(x, y) +def barplot(x=None, y=None, ylab=None, xlab=None, title=None): + if x: + plt.bar(x, y) + else: + x = list(range(len(y))) + plt.bar(x, y) if ylab: plt.ylabel(ylab) if xlab: diff --git a/dependences/__init__.py b/dependences/__init__.py new file mode 100644 index 0000000..5ebadfa --- /dev/null +++ b/dependences/__init__.py @@ -0,0 +1 @@ +from frst.dependences import pybam diff --git a/dependences/__pycache__/__init__.cpython-310.pyc b/dependences/__pycache__/__init__.cpython-310.pyc new file mode 100644 index 0000000000000000000000000000000000000000..4109ee754769b764bb611a0c696ea8364a553358 GIT binary patch literal 222 zcmYjKu?oU46iiwK5egmMojTa3xQh4%F5S9>w0YHH+Jq(!{UkTPr>m2{;N*)_5AKe8 z5AIBotQg__T34IkANfS{2o>(XE!GgCtmX2 z%e%|HJ4Mlg5nj(_xf*>8Qv?+DL>H4tCV9g&&pBr@II<*s)08EAnAyG*`j7;9TzAf= z5-l2urb5?^6%W>hWPkad;o52hPbys9cIZ54WmMG~DKcnO3AU%L$g`kg2x*ql9WL-h=UB6fE))PE*1b1DGVu$ISjdsQH+cXDNMl(n#?bOq70gh zw^$1*lM-|NG?{J*q!ksH=%u66*qFBXf%gCGdL^t^;50!xBRfTSoHCS~y{(H1F+q$CRz3@!F9fF&2Z z(7n3^v07MiC`U|Fr%Ll^nm8%wCTW#8X}>(bw)xT~jnkxUUUl=@y!?t|I<=Ent&_-! z)&D3~t_%&m)|$ z--mG0egR>|oBsTCp>Ngf>Q?O`ZG6LHn*DzJ;DH`v$kd`+IW*v5G+^=jrT7$V8Di&L37YuYKzS=rLu zT)u2gWeyx1nI6Ss@k-7q=WJ`bFd-u-d%d?fKG_dY1 zC{TOnS(#+=Nhd#(M^?+uxmhPa&w^vgCEh0rEEQ2sMwjAhJwlDiVdP$rDgmi?>ue!Y%wD!m72R^a4C1iFskymIaVlRblx7w|gA?=lnOref z&S%|W>!gzd5uKTI^7Nu$QvHYV%-WEu$8(jko1MiVPz5n~&N`(TCo?x^4IMvm;%o|) zETR*TH(Gh{P+!`Aqh5)i{*Jcj_(U&{6!(;`ecU^N63qcRH1*>qqK(wyqzm>uy~? zU)S?>{Zw6_zNrJ`O-w88riq~1d-%w)(~q4y|ICGFpL_g0$B&+=cAP$S{LJ~s&mMkn zJ?3WSst;HfK%=ahKWJNirZ~q=MRU&|In9|jmno;~G0cnVr>!G}g3xqou?*UmXVI*= z%;lU_nV00%V1@>qT*e+q4JW;Uy2cMi20#Pp6y}|?XkDs+ei((IF$W`Fuq%(UM)^W6 zZ9VPgrYnUOKo@2)ezT}RnM@#K6>|##2#!I3LSJB11oR zr;cZ`F-=pWWENbyBrNdMacc<7LwW`E!~mCRV3BS+Qt1{Qb1VeJ0(K<>TOv26{IH}Is2)-XL+Y>4mWnwHikvMk?b29_Y7fjUPJ*QA zmG!yEF`j03395+!YX}0X0-hjg#Fc6R`Y4z@r5qVV#>u;-BKTlMgck?Cz@-?Qt7K=b z>5NMrBGXUuv;a;BnQXZNffiCxx(c)f(%jX|90b}hBpN0R8zZ9;#3Gx^w1DoGL`M5_ zQ}}M7RJu$UgbtwS-2&>Y$Ybp|jH>5d_0)RX^^f(vhE0-AQ@_M~k|JEW&!0Y}d(^D^AfG*fn6~ zr^B)4)7oQXG6t)+!^otB$H>sw(dr*t7Y~$Yn=!J{JyRJMvtQ8!s?aKKW^p96XqeEr zL2@Z-yH+9Rio7aXqrz-L-4bkJU}S(JwFU5jVhJ+7oEeZ2&Y8_NSqr3Krr_cwHNHGV zh;VMa9Lb!9JQmhDa%=5VBi(C3rqD{j`P5Pj)T80B1qGCBlrs{8VKK^G!KS^cwBpfE z7oo{14#sv4gHd`@PX@3;c{w*PVtU29R?;4w0UaZc1}<7N`77iUWW}5vA`CGQ)GP@# z03{vL^>)fMRH%A8A`SU>dEGLXnSZf*)Vla!^Kp1OGe6JjJOyT%FJn2OblX*m4WO|i z0_}Ec!UE}fFGo){Ka+7LT`U-iGu4-@iwB$6=-8>pj+}#VWHvB1unl4Txu|;9V&p<= zD=qjmRUJ1FsDT6NdQVuLLf$P`f6F>zF(R$6Kr&!80Cr~W^JkM(fiGG~O!jD~o^Xmo ztoe$$YFP=uz|t}0LHa@-)x?Yh?{r&Dc7c|*PI|%+(mhkiRtgy;$MTS5shmSlh^48UR*KJ#vxA zkUDr9<1#Fkt|!?D)pylDvo1c|ysuhH*<0)^iefYstG};F*tVc4m6$--Ds~aJu9>_I z8$Xg6zh#N*w>2^;|IrCvnT)q3dRx$zy7;iwME)+8Dl&MR%Uw^M@K~sfv+BOv7^+fY z#1p~gLatc7-@5pY=7aRY6DQ6nPUi_($$vXSmIa_GaL~B97pvcbW&B8UT<1?bt zlj=@S-6*Spy# z<;;vbN&O(XI*3 zEb4_a>q3T-T~nv`tqugLX21{yB6&9KdD&LNiyX_&l{5JQY!iWgqX2nB36P_rjEH5$ zEcUe)eK8AVA!r2FGStflva)6YO+5Vy9Y&x|rZ|(c_Kfa(hP%o;fEEB~AJQYAD2e2FulwGWMZ2l#btm7|JLzkRPDA zA`fc_B>Eae*f)CMp#u-^KX^dVB7AJsiiT0DBo1hcSFR zhGvM3#G?<)cvwDe2+(_rgCD_U-xp4H(SatcbFykzDh5VZU?fH99Dzdh850#|=QxFcfTj9kf+ts$AtoO|3!xypLccICEY)@KVjFhCkQxRip?6&o~MfIOJq z?{j{C=1cEj7-Meb-4?V{^Bu!6ScupjvuAT|huI^)x7wI0l(Ls6-F!8d+3jmKkp-1! zhHqxTFT*bH(l8$m%@*~iSrdC4MVQSvxG`QF?6=U+@AFNWpz<=Tm+GAspz~!nSC~#) zCn5`fj`9eY-3P40hdq3>$tq4Q7+V0BnYeNT!&9w81hArCK&x)e1Roajregut%8WfJ z#<`*Zf{rnRIq8i`wIMKH0PRLG*|HgsObDf3aMh7>KwB&a?H0{ZbIv;R^b890q@;iY)eYOQSaCOOspO`?AIe;1Az5@^sg_h+DLk%=9# zoJz4sew~F(qB7^rqHs8dDq9$MibzV|XhEug z^deeh#EFb#>xL~@eaU?=*|8f5p9E|;aOKcaRLtdI&&yJE;ck?S`a&MQ52$;kD7I?a zLo)>;U6kqOauBqrtIW5Cda`n?X}aDg-2uUZ-PK98u~xl?op$V;t4!E}5PswH0M$Tr zPhM^L$|i-%QJQX1kq{zllU~8D)DEwz1Gt;280Cx%k?LpG9v?~?5vFn&37U7yn*&gY zl)mqZuE54-9H=ndbI#6IikD%Z6jK|z1}i6+CSjw#h1O=~IFi}>0)BNB_s-|B3W58A z^3kYom1c*=vq>dST+gVeI?D2PK)vbVYnP~X77ig8#bJ5b|7oS6P}PqV@ldjV&|Lh_ z!1F|T$$%U{8s^xQKH6wi}@ttbZVxWpI75hX=%bXp#T=3`n94OM_Cl1s6%Ne zyJTJP5)#v*bSq3seqsw+HNXSqr*dS6)U$@I0akVZnu;|5ub8TUXC2%{Z2Qi&04!@o z@kbHDjO;`<2*=L(-H3nhlpf8qh zidB$V3*og`%I6`BvAdMda1A4yDkSv0IG42`t$-;_(VZ)_6ni#bg-E85Omvs~^k6pV zv$R{re)@2$l^5u%R`C8_(?ZFSjpLA39tqU=t9J@?C|j%u3O|nZ_XH7!U#m-O{2?s?$WiS*lyo)~N_2cCyf3adf}J1BWE>k%7!F z$KZc>Is+Myg?oZK=J94S@2|;$?xJ>HkrI=oiAh0EDMLirg)HPm^+Q@q8Yw8$R!x^4h~C< zB?cz~>}m^xHc`JHRDGC@RXCl$ikZ};2FiiLI=prZ*Fp_M9uDPbgWXX!a}q8wAm=2_ z3?3z;EmZh?qH*fuEoNQDDZ3n$LxZ%$4yJ_t30WSoMn{7D4H!WYw(Xc=YsVvEhTdM_ zbTXXhGv(P-L+`e1bGlHS%M1>OjI!h$9%P!eXHS@j2MkX7jYLO;Dx126im7SPx@sX? z3URMrA&~DrrQ2Y`o}8*=N$rD4CW8c|J6UpIKtJWlMoT0w{55@vj9I|h(2`y&N#)Dc zK_P>TI_uez1J|vg*$g^kqOK?FS{4E(_+WNx?f4$U*L@Vhj8@Yo?t0!>(wB^yZfh0I z?ZK~J(VSk}D8lg^Pv%x_S)9#p7+bX(zIOXo?OTl(>9U+~4grkQhoEXF5C8psk|0yn zu&lRq>*jV~d`sWsoWL&~Mhjeyr($&@JyJJ0E$W>SSLHV}=L`#+-8+lb#a?Wt&t+Wj zpPXCXn}$ItQ`lR`Pr>aMZ;JcVqa%Cah?Lty*VD{Q&fPnWwZL9GM@LrfhPcW>I6c2u zPadM(o_+a91sGfd!K0dU7C%{h8$Kleqxib-My{NO695{$|Ka&RZO3f$RUN+n@acb5 zgYQ55`tcj5Z$Ez9a|wF*zp7o+$2$=3wtFO=L>zwp67NKOqunR*F2v#IFY#`~;q5Q+ z9>n4EFYyhC!}DL_y@*asCZHfK2~{!oL@@{29jyD1op+Fh+tKfHO25PG1)bA0HJGx zU~DdmSGRVBFurYREdJpQm_;P3w%)c_z6Bcmoe1Vzs$MhU+`NK+Z{k_i^VURk^Qtwf zsw?!FsE{VA0o9Re2H;oZ9C(lP7Ie5cC6f?-d^K5DX10K^l#8>#&K=p zz>--rmy>qe4c#4+Sp0_WoUfVXc)4@Ai|4DZ=*~I7XfJovOevq(sv!>efYqU3y_jAz zM^fdut~s@G+jYPzcR$(OE3n<^sS`|&idKi06drg1Mn2w2@O)Rz5 z+ODzP-S)<%INBFS`}!_x&iB^h<(})tC6sl z%=s~3=tqeSsQ+e<7KtU~u=@ZZQSP<3ge9&)PuQ+4b<`wBtwVAM_F4ki{}9+aeEiWG z@U8o6Asp0vqk`dK31GagJNmj%zSOiv$#P#U$acFF>+ZdPyEUw# z*U~PprCp6JMJ>;-1M39nsV{F*(Cls7G*WxH-`;)$XY+Jz&H62^1GChue!vqWF8hw{ z+99;!#ok&1wAIR+YjGQRi1!+23HvZ%OFfq2Gul#jxdY`kUelMi+IQX{-L@|E)Os-X z@4BI*bw6M0M*Vy6_U;?Hy$xgK9#xw42cD1Cx{3F440(U2yv@Fs&@h6mtr(S;BJ+yr z%l$VvRySh&%4yD56H6Ov8#rc>x?O2-OT9G%u)bRBy>4JO8P_qRmQ)F9b%~|EVqD4( zs;62X%J$eA^p2r34p$c^t7ntRlibmPkrM2P+(j)@g=1mWFDD9M8Nq(Ul|1(`VF!YT zLaDeXH#ZMw9IDW1>p5(y&G7VK#uBw&Ruv8zQFUuF>8yd zy+nfvH_^BuCu=gXUBiRXr5z3iAYr{bl4s?3;6Q{CdmJ}&1x_o%Wfcw@YCLb>+*P>a zV#5-U+;p<~)6snsEUlGcRoocGR0+$4@{E_w7q4EhG6Dn4L*(NMLlc_8^qFo_MaT*f zC~#=yi4R|9C{c0C6|X0Jsq~id(o*#&lgVdkX=`G-L}#Mfw7~{0JW(4JJJEN;W~WD2 z&b}#8x}FCzXjH;hCLIm1D}cS6QWkbS-?wbwrKJJNJPZ=IGUc6$^i0(mF)`4>5bhQv zGatp!OD3O&6=aWF@%*G4D&P+~zm6D}z^@MUR4^1g+oH0}VhEy2VR>{p1}M57`zo}R z^7xV%$>7bzfg&b}v1Bei1GD!!Mhi9z<`)@eq-8h*FgMU3OoWi4a@+`;moF_bYm^RZ zh@ko`lun{DvirOa<~7kubClu*)Ucmp>j=|A5n15FImz1-^0pj+H(JU!8`7++497wl z(xq_}aK@p*JcOD~W-*V0&lbcw@!gNF`$gy(AJEELS-%7+cg3+p}76gRb#kISzfo1ea^O)E?C*s-0zh8NB#X+j7SZm@RI=xPelKgHrK^ zQZQvu%6k}-z-V75cWhe&&h{2?;&uY?QTtFULA#*CNJ!p|`cfK^i%QytCU1ISc1d0~ zoS(Ck(h{uxdA58BG%&SG8zFPSr61A^ZE_}y{SeD^g}9XBAMEPe49Vi#(dW)nKq#=+b`X*EbVJvU)1A)pk*!dx*M9e zu6n1`YRimXL#t8St*W+y|4tgdfah^OE#D;+*|HSFsM%*~<-1?k7nfA8dr-Em3*~YI z$~J-0*mGJ@U(i%ogKjBQ{&GO?FNgF-y}pliM=n@vuIn?VNAEj71PuiA z?r)@bEe5*Z1)A$CDPF#3`Ce)jGDfy9TbLahdR%wDEK~@3p+>mqnUy0=^n^dXjlTFz zz-?0_EIQI`^btlZW&qwM?0aNnn;6s3>+xnIw4^rs-le#0L2tI83GcWx0L|CL*gUCn z)!Mf~Z-@41I0tKHtzBtIjqrl{11j@-vjs3&0%h8^Dy>Yf>8~5^vDft?jh)@ltPHF) z;`Wa6eKj?Y2NaK}4%XVs_uKag#`|FkZMPp#Frc$q+cD1VyD>+TwYKF4YN)SzUrX9M zl^$vjUNb}sO)j-#Z2goP4K-bIqBh7sq~HNcr#-aPCD=M)MU$E;EKA+BE_-O3R_m7b zf4bHwotGdJk{Rm)vRyr_vls70eJdA zyG>pXY_j(#ia)18E5G|;tjq9g370!&r&6Mnw4Jh(iB8MYVlr@yS#`YGA|S@+WB+UDg{O;#h%V*Z)IUTlOba12;T=h5BlH9=0>Yep zr&RxeYHVOfvigo>@|5hDimJjZ5(qF*{-82>Do(I5Di8sRlg>HBobw1SV2h#6oy`>r^|(?39^W7@J?*n zm)~hESv4zVq6l}qo%a&Va}352jGJ_3NX49s4as$r-VW8S=q^|_;ctVbj1WL|YvWc~ zuvg}60K~k8F?4euVRMyMJ)TTnfRY}noBmRXia8fjE>qg0~x7)cZ#R^Dk(zhSCag?2U+a&w5t}ewn;W#CyZc;CIo+7k3 zUlvPsQ!<&Z9(a`S~uHdi;#o;rJ?-Z>d&k*=sGmH)eJI+>T@1W<2(^1N47UUh=rFXVFb zb#uPNwpFpyR8NLW4I1y~9Vf0Plt{T_WAzFm4;e*5)1j63v%-fQ$B#l*Ku@6*k{O_~Wj zC5?V#Sl@}%9;C+oU+=A9ZUKx0-Wj?dOX$Yi<}1cm`rC)~Aw8-0;%&T1*?6oMZ%ywX zA@>2+7CbZk9zfiwTY$}4nn>@#Z({_$iw0SJ8R zH~P$j`UvXLy#D=0RXqIClau*kzC1aZJ|sJz?jz}de9koJC*0O~2C-^)15>SHzuWf| z4bQ4z$({*ZmC{Kk&Vt6Ff-oW6$ctBt!QIo`qE|JxsRtJg!bKWgt-K^2PH{IH~Vn%+u&C z6r8{js|lCHIDh2Q*e-2~)E+)#yfel>sgpuOkR4W-;93{!oLs6SCJTC;!#^B0C2P+k ztYBK$xeKspys86i2Do7ZcCk!hfmbh~XGjBN_MD;7X>g2|Tfy)`U4wu|jfr4Tq*9PT z&C+Fcu4>GxJ`4j9jUWWp;CzU~FRD|f!P!a+;=df!9^ON%oDV1|WX=CoL_UCZ*b;UJ zmW&Rq`Yd!LXxrFL!1@Ii?fan@8%pmj>({kQV&_-68*{oF%Gf(+ZdFz+uLiB<*R@IW zI`pecZFZYz|83EdfC+sOn4Z5iOn!d%ICS@<3B6`<#@E%T=PaQ9ejQX^)&zCH9thi{ z(y+$WsByj+PsNq>^n#jhoHa6I;8Jk~aY#7~+In&_oSU31&H)DWxu|!aFPG*|(IJ<^ zUhJe$PP#mfEi4~Wjhw%(&%dEN-^r5ng7A1I?tYq-y--OID$h)WUZbR^{*L$~-@7xKyPx9oEXH$@;-^;TZErC*wd=-1E z_q$Ew?YpGhM{B0@?Zr-%{vN-KK)IWx zk#dg@1+BiVeb$qzg254`1-BYP)*VJ0wYapzOT?SOp~rKObM_;WlB1~~o_PeTEYR1?^{Ja`+CbGQHA9`|v2hZ=3=O=<^v#(1cl(SySk0a$N zNjc}G{0vggOUeZ=oC(%8aDQ z96`zuQf4JZ#ye6*k&>4bHChd43@Mi+MP@ScoK|xgu~T?=neQ;SI_!evlKG2wFW5fb z=kabX$W@eFGTV{M4RV!`s}$s#ms~QZk?R}PtVZk`kn6=Dmm|4k&f?v_RI?Vb??5g$ z$W@kHG7ItUBY0Pl6yb76`C+77krd&ENcjn*EHDMXSM9~+7_TaQ6@uwFjNEcLX-&*I zcwd8XEVp5My}g{k)^XBa`fz+HS?p)7+4JMVDJEg`x)$HkYcc#!5#)9LBQps z0Kh8W@;*e$Ynb3wJ2WCrH7a4mI)t;?xdJpJFUyFv+Zs(dcM%BjRNvJIM6Q%*1jD42 zfK-}*8YC!b-)6$w6TpK&wI@r{>WV{v+gx5KOIM!423Lav@`SzIovY}zI!HmPm(xiw+lGP%|Qy>}I8 z-ue(}EF15S5X7BP(&+J4rPd>U6kL^##p)A`qDK6;O|1s-;lbS%;J;qNMX?C`s6H9A zgo91&I3VNTu0en5tL0QU_j$USwv*^JnHQxj*Ig6NHjYd3o=E>I&|EAwfac6scQZzK z%bM`R+!M}ygb$Ag>;Q|yKJL&cOnc`EWOG&SxwZ0hJfVVfcewOc+9~f~{;Kfcj(Vvs z!c}jvSf1p%SGX@OOY#)?G3N>3t0(1VN%8itchKM^SBt>pM5~&_vozrOAI7;Z6D9o< zh}5uFfpTrtOs=B1H|F%#;OB)s2c#Tgh?UJU^2u(7X~*PAsoM}4K{?$iB5t;6fC9hf zW$bjYqzEGHyV?m!L1{^W4_DbxWo7pi9$R=8yz<;h(G}AU!7h$mCv34ca7ex$qm?gK zp2Nb71}@LB9B0!aE3|gHpHEiiKgabMeRS(i0-o&y) z$8=5tHT*)0j|bZ6Vm)bHcV1wSK~UF^N~3OeGldfSCXH?mEy;ELaffza=Y0$&8BoIC z6j?}~f0NVdCPnQ{3xSCAn=aSzs>p;UZ9wtPG zSOHCO0Wa%eH9>&oVq@ts>*3YV?p|Ft=oewW2(-r1W7g=tFq_0Cpy91*X+#I0ybD^| zSQ^mwuYz{(YG`t7$Ilpr)=(PI8kL$n& zty@EBKx-r<>$Y_@G+w_PV2wg+C_QHFYaB9!_KwxifK_EmAW_FWIMI>n z*GLr%0R{P20^(&oKAdKNR6+lb69hR)>|Fsrr2esj*kqBCba8CYXcTFaya7jP?9oW= zg}!@6fi%3ZfPjQvq8gy?7l$l7TaP?Ko;%DqZe``}mV1%uJFStaQG6y+;_b;K?jTZk z85%q_h-)`|fH;FzeYK%+3YOfY`t{=pn(KF#B3;;A!G>y1RF|66!aGz8O>Z#G?Y5zd zey`xo5Sv8`q@aqbS~sCUKme9piwF~`>WBZq(ce04Qq}i2)Q%%QqzO*;*u{+3`-llm zCP0(M)L=}24)j%I&PJLX^#%sOrzTPt$2jVQ6;G_7&fg>cRYZ}Fjqn@bo>*DCN&|Jg z<=d1$W7pPA- zmlXe+m`G(dVZqQ)x#0Ff@{VxDS(=`9b6iOU(-Tb#3SyrOJ4#+2+7)8XiTC08iT^Dnie2UpRq z8M0IA(lcC7h5h9@9FdUS;OM?Dt_q?fx$L#0#h3UelHlUQ-*n|N|ESckl~vP0H)k!J$NQft+wkv@uCD|>JM7^YH_v4EQp%R`McrTAe_h|ZH%UN=KRf2v5s zen&`zQzg|{Tn%4LY8BxWuKsL5#{hV!Vsv;q#u3PZ!LdPvLxXv-T^;3M7|b4X&z~Ap zX6vI66+y0{V`;K7Dj!T7W>JJnW5F752snn1!RAUk7;_FCOC3I}%u31#4rK=Kh0aHY zj}9LrKA~JtP3mK*4F57ivln->Y;>#Z6_hcE zd-ibqdk=o<`qeN^_VMMWfzyB_HxKk`)py*cB`CbVK;5Mf7f&>Yb*py-T3)?PL`1z! zSa7{Fq7v8R%D9y3g2!I2TyGQ6h`-!5kBO5LbJ+N%CWoa%JkEarGUq=s_)iRehPch# z)m%2!DQ*k&^W*uh@eVq{i8eakg*g{&wFFO6ukob!4!t$r>F0y?;rtN6{|EzG{BV^5 za7^-sFGbZK<{R-S>uEHceKSpPbSv{+xJ;VfjGoB13vn9bRM;zVQ{QILp@L~WaINT- zCpellNr|TR1^Jk&-t@pdjqp`y#%TIe=D*$^ngUnyxYe2 z3EeW{20r2PG7LA<;4Yqer%WHU#HyRTGS2xh zGm--0xT+?~F_vqS%dExDT#wT)68D@W=W=k@T%@a-b67IDIHp*xcbw&)ZhSW56rKNu z>AlJEctO|CK1aL5=WIBcZVBFI>;wD_G!a%pF$6D9v~Zf(QA1V(W$oG6XAt(a-hY0`KFj-D~?<;1;ajArw$w-ZJ$ zjpwl~m`A&qq7TE#sSCMx@fRUoMR(B!6b_;=tjGEh-ldyxCWZ55tSgr6Nf_G8@uzxN z+DPiv?uN17z}~?wa6-$wZAXpqnm5KHV=NqA{DCPnT ziiT_!bgb6!Glc)|7_14OTNb4Net*z{$i78uJc8oiuZxG987~HS)R1F4s1>$GAh2co zY0ZXke}f}hbwx-hY;v@QAMR|e37mM?zQS#f)@UY4q90?hCK_&qECgsmI#Ks- zjX+@k___$VJF*xc2#2iLOQo)XTYczu4C)Y*8t(2tzP1q38 z>swWiL8%C}VlHh-gh$k?}B#{3$aDA4zx_n$*tTjcJMz+@iWMS{J zz?#^C-X(RFr@AL}y{>is1Yn!u8^w!nBEC-{;t5?gj*G`t_}Y+%g8fZyt(PGR1s{CM zJx_O2-1OSc$ia)qUcQd(P^2kM*A|w9`A<<38Cu)?ZsiKUx!m$A<9rSUQ$Erch(ja{ zkvjY{&Q04*zkqLHyFZVJ*KXp11RGhi-2p0-U zDACr(cn>iC^8lk7M8w!stgS-D1Sgbj>!YG8#z>hF%Vo@_$NP%y5KjE=j@0kW{u zno37bD)KXzya&82;Aq}$O~wr<`c#QUB5yz@UvhjBap%_=++ZL}8wsA|V?TnA1eFN8 zv$>39`=k~Fh}cl}UnOsik+IhrdCi;*=R$LxRo?)Q;vMMN=OUne$Xf6kY3%~6kS8(8%Rl(eru$J zdn1mtd27?JiKyEndx|m`?TQx9Y{%uCqH_%S!i>RAV|n|E*8ZeY&RtdZ<>>am94@ha@P~Fr9RnaoUz2cg#U#xImxK2dpDM2H?$?( zwWQK*q|*16w|LkF{BI94eYp+1)8L(sxDg0GF%)9r>%?Aq?P`y9BSw-+c`y zJaJgD2j0tOEk2EdjCgM20ltvVZdZ+R1-{6==#;uX{wDc>=6p8;N(qP04(E?D4snQ! z45_~B>)dUwUI_M(8TmWfN3P}?p#B|g7~`mRbhmgcmD-@zRaZz3nIWPP{5c>0GlOs+ zDaJ=k4$b-R3_i%<4Fud<%9q88o@)0tm~a(%rBJTN)rRhPf;+%mkgFIsR>za!W-qR0 z@Eki<&i1{4|0A4`X+%teNAb77=A1xy7s6hB0O6Q^zup78bprP4IMPhL=L=?kwJ+cd zvWv)#gH?>q1X!8DRFp%3fsv?3nA=-=M?*=LqSSWzeF4D-Fvli#KOf`8?@PS*2y+AT zhX-)w;&-{`9chrmG0I@fxxrct_p{+7vmILoH^5DM<;pM5==UQUEz*S84`yFm@&gsPNT04)JXC z)p_TUUn0TTM&p(DMRg@1<0?V#mjdH{jS$&oO_ie&y!b_*xo$D4;VLo$xPi``vK(;i zAz81ZU(ic9W!`RMw#$>Tm4^ok|I+Q8qWr909t{22g_KBkC?-fLN1-vt6wj5?1pxWS;Y2jyhXub~E_ zqutS9ec;?%yM13Ju(uHa#eIxF%AgsF`|;v8*hn4Z_AF;qBE%D56#{Y$nOz*bs*RWf zo@~U(#~4Q{=Kx}o4jB=7YV)6vD=yOU0XtUYY9bz#rc_E`KJBwjd4qE=0W{pC8+q#V zpHQMJDBm5#dN`7(SN|1Tf06+N5B@B|VK7N9hm?v;PVse}!E8oXnif9}y1MU0gnMwt z+t?9)+cb=?Vu$#viJ0-X2EGZ?_!|S?zBc142EP4q<1K6=zm;e={@TDdkud(sz}M_B z{?fqLOd4M{@HIP)|7W~ilci!IDrRcEVU8;3o97_2|#y?dO17R+9aCwO2^$H@TUm1DO?RyAbhlOJhHV z1Rx~HwGsXa-hSRk*riuDHXvN3mg>kJVB~BiMjBLc?2HN)3=21?PGQxDI}GX!eh)!CflH+3a}HJm4>FA`)47|$Lk#X^ zu%7|9lH|;b992;VqPeP;^BfKh|IPpdn!_Ec6mmjTu~m^7GSF#yl=E3?2+D+im}$Zr zMF7x5qjY!cZz1Y zoQwV7{sMvoyk?U3XfN|$;>iR&SBzK^-XmY<7rXCY7Yvr)hOgNKPH$scmmw|6qA?90 zG+%oZ-#0l^!awK&rffJLN79)$Vvb1~;pz#^iLXU1X8jTdrHH{UO^U9G%jp0Zz%LTM zlUxfXUlX3w$sBx2ZFSli&m#^`Ex4InL(F-F!F2{#8K4NRA%o6@qo-=`Ak?SDP3vcv z>sJtr$ME+Hl;Rf0A05PB>V)eT{IgtjHCxAttJ&PVTusCFDd!g%Az{F=QcoOG6YLRG zPIKQvU}}j?zI&58^SH&}U0HOc!dHZh3E^Twu$T}kCIpHJVPblMghGS@#OyNuui0(9 zY4#XjGB+51VfGq-Zf-RH%GVe7$X*lEOYvw-VSIzzLUoj63 Q-O~ei67{e5<$nGD14;+Vs{jB1 literal 0 HcmV?d00001 diff --git a/dependences/__pycache__/pybam.cpython-38.pyc b/dependences/__pycache__/pybam.cpython-38.pyc new file mode 100644 index 0000000000000000000000000000000000000000..e0d7aa7c0c3f4ee3274e6242ea1c2d27496d234c GIT binary patch literal 39196 zcmd6Q3v`^vb>9E~`@muWf*?roC5ie8Nd#5|mrqeNOw!_0q%2Z==m`p5Ep`^b5{q5% z-vvpm7I6|P74vYarjC!Bx+$0@sg>inr;oUKoa56xnx<*pI&GRJ8#m2KQlBtAZJgL~ z62-CmeRuxX0t+l)IXMM|AM?*YbMKuyGxy%Pb7%D8x^+ndKR>_wk4Ap#*9_zL`4av) zgqNrBTY1(r3}uWM%2Za-95LmuHDbx%$Vf!~Mn|IZH#QQJzwwbc{zi%&(;XunChJCv ziRsRfPNZWht~yjgb*iN5Qr&9ZClVt`)g%AbtM#f^{xP>pZ5ZiReIx7CMzu*jpf;;5 zYOC6&wyOu#L+W9*Lp`EwwNv%0N7ZBMakWbgs6mxdPpI8$NbON+wO8#^`_%z;P#scF zs>A9jbws^S4XgL7qw1JCu0EiiRwvX+bxNI9XVlqG931IU&y1|kuUF^P`GsEftoq=a zMt(!SXQA(&amPH3-+Sg3L!H}VJY|fmf5uSHsS9r!>cSmsWFy`#s^{_ce16lN$jAeD zzoahX{j$7oR@PM`btU^(mSH3>&5dQJ?Qy3(ZRg6<(}k)%mOXTM-}rvKmTu;qYF^pn zWyiiWSDh@E4%%aL_PK1Sl&{#&%@$O)n4g(RCr`~fPQFwv&e`P>9%`tQ$ zlh2k*_N8obI`7zf?31$-_IRPB@UqL!s+)ye9+Z29U4THEypWxoo~_z17H3sD_c!IBW=@Pmp1p=@?`0`atUoLSF2^O)h|qfPr0%}%d!>w zmQ$!!^Cj||#qFT{ke#(JO`f(Z)ma6m+vRa->7{aIrcfyq3sd%U*)cm?mD~&2+xDem zwv?N)&zCCILKVzmjq}sfv!$^@xmcc<0}sy36ejYee6^6P4B2O$JecTAq>~rs1e5MR z1+Ib`a`kk6wpz(eVi4$x7(ADp@`RI}p0)>1pE+|W1tLr61e6U^aN*LGt6s@;^1L)* z&OV_EGv!L&KIxLfnx_g*E~^UJl6_?oBX8O+lxz&CJm3~F5N9VQZ4A{)K0jR{7SvYR zxrz;tGkH{>E?0yx)1bG2A>(obV{HZuaO~kdspNf1nTsj25z9JGb`FX7Y-uK&o2pqy zbNC&5*E)IBO4fATlEHupX4RicHM3iwK)5^FX z#iQ1H?8K=H&s@3s(hD!YdhUa#PhPBbUO0dH;?;APj=fTQ+Q z^y-NVoOILKYP$CO_K9LqI5{>~1;?wCXv}nWDsRuuNcK_)gZ)lEtNK$zNq2ZI^9SPr zpn-1;v(1^aug`*O7<=F>haf(%15c8eLNTAVU##TEXNzrs-kQXy&4GX_r+(Hhl7;G5(L?-XdjMY zflEC!J)4`f$Fmg*0hvisrVVgH$mXiE&`1FnrK`YOAkE*-PD3LNK_y}AvoSKSQaKw13Z7wu0HM+YDP}iD0tX({0)yJZ#c0nwu57 zn6HRxDcSo)m;}ZjrhWguevZU8!23&ODB)_hUq&!zE{NT4-y*ZJ#R@*sBr8ypw=^1fJtfu-dVTQ>yp z=T0aYmKx<$en#}@qHpcgF}eWyzkmkL*%O7Eoc|QC9P^^X(EWPs2P#_2{B&(CN(pp7 zk*=MzuRYm%s9ng;%ns z0`dp5TO@QoYzws+fUuLww>&24p6K_@t^RcFH|!I3v7imZ>MR5{tZYab`+zqI1SLPT zQ)KLtz(C;4gtqYvUfq&1fI-q?rh=Wt0?5Kdf=H{ho1+3PZJ%{@6qIhZn42wTwHPmn zO3W1|vQ9=bD_uLd28%@3MF0P`RPuGM*UP??8-T7P&Or#bY-ui)7WtxQqVAnGz2So% zLM2FRJ#=%h=*MU@xU+^#!$;~5GF1!eHjtO-U@QZKEi?r z7fVvNjIf2{lJPzba)Me zjoqQGlyvQn?Q2I`_eWdRw#p8n{zG%HI{KzGTpPa8NbohBp6u)fJ>*&`xsg>XIj?=8 z8TS93u9KOyCT-V-pEPQR>=r6-t~@J)vbCynHF&yqU=8C?I%!|J_LzO`eXYmf3(uXo zs3p|ev|6>8FQseWkJanG)_|{`x$H{5!P6H{z+rJ?HeaH*NzXiA+_qivwNh>;PCfI& z1$Y-+=YnXUW0<3R*T5T4l&;O80sC7w;PQ(np7UC8!?g%1Ijl6bhGwIvvD)-D$Mc!X ztyb=9v0R8ho!!9YbZbq1P*ds}52DE-fUCF4dy08@(;WK+tnAO9rbcQ`*UK*gD0-_f_@=TCtI4x+k5sOc#2laZv1U$t5R?GpIQkT zny`=14m|l}J5HQ$Qj@Mq152=Q0G)ZK024_10!AHFY!=J4vUXOXH?#@dh8fxu&|m21 z8odA)0JQ2dOgOOr&{Kzw96WqT^CEa{*9xkGYIdTM5gSvMJ>V@?c636yt3@(97*sCz zy_NrJcuqqJ%@D_c%OAL@;6eLIOyL8;6c)2%)V?B%RISx9n0%ieXw4G_Og!4!x5~3c zICq&fi;iu>OpqP2Y6itR2+=~PJl79dy8@Rc=FY(bW3XeS*%%a{KD9RSElh>U`b#3# z;Hz|J+hQ4N56V>KyrOUaHi@7!rqWipUMl7x!?=j&q62P1Y#GpWe=>d0lhod<7S3ac zV%`)gZ5XBJGKQLO!mG6jXS6k)-R&6?p%sVg32x;+uEI53p>H@CU2VvuI}XnE65Kr; zpBMuU9R2h)d!kR~n0y>7mUB~?N}=ZZVcLK$R4e)7c-lS_TH11+guv`RWFI@`;-fE7 z3tPe12DnV7#qRU9S*QUW=~_vK7Zmygrja|gwD)v~&?>bXH|(^QUY7S1N`^1J|EkcT=Ccj(VE+w5G(Ze=Ccv z1q|(}Qu$UJN;$fxvE~!SuGi7-SfVvUHwrvC6h!PU4<_WQY_hC!fh*pZK`o0N7)98^ z!-BKPU7sJvOs+y!ti`!s%!onrTZ>qFF>=CoW=kas#vCL)Eh2ZI z2cyY1Nh`MAN2;{nNWv~_j;1_!{4iC88D)GD6WQV z_!Bw*F(0udh3#m#s%cWUmC4h2xR`R(|J=urL0c?f+Z(jcmc$E6pIf#_q;oP`D(v1i zw90g8qEjR(C{5SyU_%ufq&fqRYMK7tbk&VPR}}OfpSZGhiSl?Gqya=?PA?KCdP`4l z0i|nyuzc5RQBM#$48W&+)gEP75ZWHDh`wdpvJP}U_Y`xJv!y9`$HWH&uEV7W#~56k zt7u$ynuC-bE8v%sxOb)iD;=BCdgtABYPq|a?(Ua6z&#@i#O1J?k8!xj57bVX)8Y+P z(R!jnuA^1>FN?4N9pZpi18Gop55;zy=1%3Yd)|i`9Cq`j+W?l;ki-cH zmkQ|r0^%Vs5Wx%bT?mDtR+7LF&erL62LB1D$TA8#u087U4MQP^O)<>*pgFWA!79iJ zmCl~?EIJ*zQpSK`LkcAXC-^-kwD=}D3_(58WL+`VZa86q&X^d_!KWdMod>3Jaaw@R z!*)U;OEX=Z?quwY1o^aKrGP1qF{~!=Kf0bojYwvLOw9_nTOl6MN9jrx+oD7578hWI zUXyzJ4Y$gU?EeM4a#^6qTkS28dnc}*I)5I!V~*JEp*imgK!2K&(O4GCl2U(b(_9Qb zN0bZf)@r$26QzR^i5Tp}SFop^VX3Fb0r zp^;#wq@*N<_W9Gg5~lWaURTn&Qi+>k4KKB?(ZTA-(^=g6(GVroODy3nj7Cq8!^N|P zE+zD*z>)0H0MBxxhe$%7l7ylivW+0UPeYS@O%9!QnjA@Wuo|pJO*9kI4mz83DM_tk zDSvt2Snt2;?sMT``QH$DkTMz`6Vkk3si6rEy$i6AiQibymW4B<6U4UFtK3!;wWCUs^cTrY8AVCIfN9p z#DhXZ9P2Bh8M-#$4!waPNwLPjsE=K5$;)91HuLwb@t0;2eDs+lUa1eXv0i4j2uiG{|B?%=)mtO{8paE zW5TE#qeot|=FNGlZYpCcQF#b|&8Y-j;a15~5xhmW7z+}deb3rr)bXo2witiUx)DXt za>6+ZFisyHHI+Q}H~W(Wn`)MAziZm}w*lk3<{sw^{xW2=NT*6F(y-F|8c|M*Mpr0U z`5nW#$O>2YPGXU{7u){R*$U)OzEa%_Z&3l!H7a{>Z?Q0j;9q<#9ZVnDi(sbw9>$+$ zC-RlO<5(%~Re46Ua>J;iqlD8lbB*Lt>PU9y@mV0LFw!)Rf6gT&a`j%^M7vN4WgBim2$DCIaRW2>+&mfO!PN<8Pb+@%Y=3PcSV0P2-L^(uwps)g$R7 z(g=>1bQjVaRG*}~kw$pDq}L&hFnLM$AdPT&Nv}s5Ve^viMH=DrlHPzc!ssR4hcv?J zCA|@8gw;!W6VeE;m-GWjBg|gXn~_Giy`;Awjj(%3Z$%p6_mbX*G{W#Dy&Y+U<4gKM zq!E@c>4%U;c)p|`MjB!IlHP$d!u2Kn2+|1Km$Z#E!uKV;6KRC;OS&Ivg!4=KQKS*p zFX_jSMtHxZA4eKt{*vB>^fT(5qz8~bub!3kAkrUH&q+Fk^aXWM(oZ1$yt*Xm-AG?n zS0p`z^i}nOr1v2GqIyZvX{29PuSj|?(yyu!N$*4YntDyr`;i`1ACmL|q%-PuNgqTy ztHvaKNad8mxOh_KF%u^mQL5?eFovHG_&UDQJVGN5@S9jx(%b(HoHkmNtEZ)E*wc}shczk7>f%*N+GPv@;cJt2Y^{%%wsw(5Ub9M;fbk~GqL+1B*Q}dY z!Gmv!c3u^Etu%CN0{_Nzw5D@)db%}f&Gzadqb97WiK<_Bq@DpdodVB^UShi|N_Ptb z4oL2N>HQ6?-71AKv9z2@mI>{K(tDJ!q;KM^N!4E6*%mRqx$rJxp9BujLTS9LohxT^ z)L!1tMkHu2P-j)Xh^I2|S{Yakc*07+(<(3>sawjL>U4g!8m(J*Ow*{w__rFL>a09n z?Om44Kf^Bif{7$npQ@2hVqu1x7_2@!U#qODvo05v(Go4S?qt$q| zYoVJ5t)@DhUjdYkYG*wv5E5Gqq=5m@IyJN#A6m9V(kpS-aK2WJ-37#IC#;>X+%;4e zTJSa1eb)evuK~h3p)7$^4^wX$Z<%jd;}OxdiTPMPc86_Wr#8&T(Z)F1*f(XoZ8}Hl z@oLXq>pJQ`Wz5IlHg9ZI=6r|hn@`l^^^SVtyR2=7+Bg+KSq12uPOfC{yu+5n7)FtYjQbI7-6N33Oz})HK_Obs!00H!FXo4TE{gsxrKi`su zWVNrJG_^|4@UDt_sq8=&QN_FsP4hnhwqtcD@MQ$U7J+D z9l3$;#EVf8ao((MRgVxFMuoivqwIQUt}uUL(>;#84H%E6*kg(L_4V}}gUH>ct+M&v zx&>IDulL@yFh8ujm?QJL1}(M3d|#JR9*dLQccCdq?#lfi! z5krZg*Pn8d5X*TvYXHw1_M2`NxKj+*BEmSzr9JuS8H6g+FizXAV$XAe2LrRV*z2-f za41}ATarn4;jUM-SWN<*-CyfLNQGWq3i{f0y7aln$^BAUY{>=&LMz6GW#hEy=3x5xHtAxX5CQj@S5BobRLZu zm%y)&`pjY|xbAB0+{X|ENkM&dIR+@Y9-C`)T<~C%_+t_B#epKehT&vBJ%OE*6}%B_ z7W^tQ%*e}N28b_&6CvQJ97ICU(HrwD8s>u@A|Rh!NGCyz1P82ud0Dj59BoViHMU@| ze~4+JnJn-T!R2lYx?5Dh8?NPfsJT<1_acNdMz(meCG??d0)5#p>vH!CfsV@8; z#Bb$KVLMcfs(BrH>y9;VR3mpy=-Vjt?f2DstHxa%9H(5$#9lT^kC_JlRpYP_qL6X1 z>m6YW;Z>t#y<*fhSIq^;5J6|u5Z3;KzH@tHSbw!1KzjR4b}G%^i*bBX&1k% zjdk*26ZDlKIxFUvd*CB6mb*)xwGq^8yS{xv+SrO>w9RtH{I;!F7W(->>(k)AN zA#mMX%=7Awg-2*{$T-=yU}M%8=yl?4>W4bvqHh+DH?bkUw94Ly-1s+u8q;GeJl52D zLa$>U;A=wdkkK8*xP}dnFB@P##ndD7ab?3!wqf_~yxtExH;U1@pi9*|w$eg3K>6RP zN9!Hhc4~&_QxB*t@6H&&WDV4bZPE6cSvTK~RuzRP)?*8g*FmrTUQeoB+AdWCcPz0=lk*)f zV{K3x*Q23sN?Fi@8iN`l;B=|M`EJ421+SUVs`1RPt9PrxtwwzvI1qDwv)(23Q4{q+ z7p+p*&Rz8`Tyc@&+k8(;IO{_#VSBLdkkY+?*eeja>+8uUq`T`qYB$D8%sD5pdR-ME~t@Wo3uc zaOW*jBGY$O2Nnh~yY@mSZK@jeP1JdRw|D-5`U4B8x~xrJHl8#t8&`~1LYnaH4(BZA z@?z=A`R>bCJ=^!yAzPQqS1Gwx3EScP(kkb4$BmCJ1pZ^}iXhbijkRTxO&yh1Rq@-c$P zNR%-#sfcr}DZ3G6&`hm6yv0>ZczawaD*#a2(!3AmZ-hDP0kL3U9N#}c*j%aAPA8Ku zz(fm7PH(wHkL>=xQz!DBYFZNj;(#9t*xOmuIoW%g|elH)~)w!R*qwpu+ z_{KNxck}h@OnnWH5$ois`&`4_-;2kW@oO|<)7gnauGaAe`D6E|o@hkb=?x2k2OTPZ zvos6VMtKMVkq=d;5zDYY8~VL_h<=Q5$i`LR@PM5R!a7BYABaV8f<8VL-N zBJ-TnEZl)`6sdX=Z_aUMMNr#0!Bi|)L>Me$5N5~R$jCVpEZqZFM&ph=_S`KvCeyvYQ;U(Dxc8qt|D+g8ICQX?5y%5=)lI8KCA z#CcYtq#DudWyByP5F<4sE$f(`QhaH|c*jgWwYkx8`pip+OMcFoV7(~!!5a~**3jte ztqfmBdc37xCYQyL{E>|kzJRk1`SEI|TF$^W#}pi^j3nj2{tQr7oi_n)q=Qkt8GPx| z`~L#3;k_`DywQnh4kybHgOM&UL{_tzf*R@eU-)UA|0V<`g8+(=M9#@+I=e93v6|2? zNf{r~=_A(k_&bqCcfcw{NcBi_8KExQbIl;dBWPFS0weP%Zx^%(sq+oAvt zD2@o=y#MjE&q!w+JQTU>yo6M3U6XjNVO!ga3l5I=U;&@?W8uH7qJZv&!WXuQky}(T45X5x9%!5+SgrM&kgH zcM^ps`XyT3ryfBMoMw`$%s zu8R*}m%Wj!Mi_fghv$-|#p`O(`+e8QMDH3{rNmTJ?E6?4DxknJ1neJJ0ru-y=WUI` zewj)-a5nKDJOp2GEk_@xhDFT4i9HO`q4Th6HE2iPNQk`?)80}&pExaKN+2}iCel!7 z42?@>g2~6(;P4hbiJfSyyIL*JoM#X=$F;acVTW|P!dav{6c~kk!<>1?bbgRE8O`8| zN8IO>kz?U+V1iX_d@vv~qLyX$8k>zcgb9V)W{5P3wXJwtRzUcma)?dvf9TJlJwRn& z`3lSf_NKDtBlwNtH-_K1?uG9**lSa9=YWc$2h0Ur^>@!a3QvvXo>Or*7nV43JFwRh zcW%P=wZhl_;le*&_J=qIeZAjA57$uY0V(}GQX1=e$?KB5A4l(O!HNa_0!Ub^%eM{k z?HBNEoBwUQ`z^t5mH(n1lrn#WG7tG>9(K#@KnvIL``cy(5fFaL)+y&eojW+arG!@- ze|H8@RR4sm{nQR{X`@^Fta{X~{he6dn$E@2R)O*%udYCQj8KsCxDRa?t}=Xh0=xtt zc7WrXTzDT-11`MZtH%KEhf9M3_tU67$l5Oa&*D8L@1*t#7yj;L;Kvo#Bhv1#sv${P z&e=Nn;C#Jw1b8f`2T-K8;Io_^c;DmWvz#>E(>}hvE%78>$O>im`6c$d-0yY2ekZ;j z@V_3^Uke8O1c$Ijd=jIG>(9eS8Ozk6o^t;lIcmIRU~h!6AqnieF{(tacC2^YiE)<* z*G2OExc~mJzVgw3@0RzI{=0h-r2pP87ehYdzjw;}Isbi^z(4Q5AC&iJ{rBzi{z3n} zPu`#N-_b5an0W7KDR=(9tpi>L3*6{b7ll6&=VSFw^}OWtxH+Fh&LzpY?B@I|a;`|u zRX689BIgCkdC|@J9pt&t`Ba(B?&Dn&U*Cc1u&Dn*V4@u5? zH|Gd)GLrMUoAWGkvXV3A=DdcSoaD%igg>T)93?q=zD1pHM^0XH^c;#p%6F=9$ij5jCL~9X+o{mzf@65 zX}J({2HjFI$404Xzf?&|=~*9hF1oWmhFR39%6_RCDW&If%&EF_Ip!2m>V{v+ky3g_ z#^C+xR29h)c?j>!$B|Q&9FeEU`3dCAG6#Qes#^;Ybz9B73GI6t{&2abv?d-Ne6K^F zR%6&-@2DoQcbrsjd?G%dtajl`@{Wa>l2r9=+!LI?=$wTDIIH*V5^{_bW{N1X+D3VW z1rj<)s(1M|p^n2{0f!NL&VzjH;KSx)J0Fkp@faVw_(<`wn~yzw?Bjz=OgI5xc*zeD zF~A)WwS)U4mPTjnvybA8d%g&B$t|+a-fizsIS;Y^2b$~44I0f5n6(&&&J!pDB;^5u z+T(+Q>Q$NYxQ^Ecw1bC%k;p^Q;?7~zdBNGq>iw+#{w38VwjR}+!R@vO4?RI(K}>nd z877z^f;qV?m{wGPuBFggkb)qs5Fwf{RNXT;xJS1O1yXzW9!RB}BP1nFQjRSlg<$Me z0;*kUTV0Muv}^)1TM$=!by=9J=FW1U_bvg=BMG6#vhm(9K-$?)9-V7fYb9QX!POWq ztbhFR(u{x2+-d?JoSa?+{=G7;3Cr_-BRH(Yc>@j(JTc%+XuWa@W;743(_s?6IrEy- zO4yLh;qQLuSDPDa*^V=dpJNgROQ-Q zyD-gDPdIRiJ8Y$i@(u2;S)TvWAN9Sq`YV>i8LnQ&`3N&EnQdkMiF&|cfl&} zUDQ4i<0$;#D1};xo$8*loX`3=*6HuC?Sx9$cEBgL9m><@gGTDP===y;==>8tev}V# ztk&s|tC}@r*SgXDn%#r_M@5Czy84Ga^>yuC6E2F;lq}h5k!$-#Q`lf4the|m^1)?X zD&~~&rx6prH{$Sr;Z}l1g#NrpbsFn3y0YBha2#d~Ezmj`LY|qO`9m(xI7SWe^>kw) z)(h{>3wZV-Ct=nO)9>BFZ^O`Xj(a`q0)LJR*y#pLz50BesVp81^Q1KC{yMzer(eyX zv$$cNbGVdqKEy``kNZ^6_es`$D#`mcAEJywuHih1N9`G(qfq!E4G$4?Ob=V+U)$%3 zB>iXDKGgJ%m(fLl^!#>9g#a~I9JY@v0pvWqlq}c-LAemL=Gw#d{sTcVNsYq3TGGgn z_HJrqb8SF7xCGiGOQFdT8?Ru9w&vP^*4%Wk%%K})5#QjO@nAL6)>IqN4hO{~HR`;7 zDYU?}K$$RYO|=27nUSPz%Tj1Oq3B}`Lu;x%Y#(SIErj;IrO@P5od+%it+_U!HMJZS z77w*Rj5jy3xgr2Hkx4+#;8L>W=wd6fnraW@o*l1{`*+lNVktD9dGyPL8rfVM(3&~~ zH5NBgD}5JK7*JEi;YBuP3qaZqSqjLv2B?SZ0FYmPF-6)2Tpp0F7y$CiF9sy$l?x|4 zNIi;l#Q>0BelZ{(81>+^0;DT?{}%zjB$@q$F9m5?EMhhVqhwth-m^c9v_)B8_cZrt zhmco}6?C2EeCAQ`d$$>O=~UF5=GLApRvxk&X??8{i&YOxk4u!hfD!OPjbIhy!3@!g*}Jb9W;3<8<89z%J3hK_kL(_FDHu3f z7+8Y#kk`RS1nCpHVg7jah{LFC;@+Y*YSpx~7IfR<%p>7;+Raeu{ae}t`oz`1rB}Jz zeCEoP=dTPS97Z90U2hTdZYyld$_h(%xzR0p&zZLg>G+vpdtivs9D`n$H6vzoTr^Dj zKi+*r6|DVnk-T?l-Fe*B#*3mWS$%b(jZ0M|U>neFHOp()5pjJpUo6k?|Fe=FSi-wj z$j<8<=5U!9_J^l&;6ZkC#}S+2?c3t^Cx(T~UO8HPo&PT*1ncl`dUlHc_opGdps!RJ zFWg3T9u$}nbwPUP7~wyhT8^4#{%CeB?!AMW4icHcwi|&kPFx zX(A1I9U%!$l~i+eJ@PQAb%Im4_F12gKJdWA=yYw26VL?%!vlB@4iv;8b&^V8Aa|;A z_56VLFrS2~@JkJzN>h~4_+aJ`tKz9G7A*e;fn(@YDuwNT@E#=xPo<6>(+(r;uZA&$ z@51L3Lnnt$5ub3*CzJMAD$D;@L%0bPS>Td3s(`eUsA(E8m&-BgH8&SI+Rx_axaH_f z50!Isxoiz`2SctoZEAUZqa6}>-NWC)wT}rOsXX<Djp)04wGk8Li2upf45m(A99|B!MmibLCidt^H)fr`L*l;)T_Z{N8~SUc%PR-_!+D3`Kg9=q zdbsWYI5ND?P4o39_(h__dYYYM-%BUk({x=i2EzT9uo7mE`@0K&yYU-5>3l0Gs5sK{ja-+uD(|P62u%y8VK@wZ@<8tWboQ5 z{wB8SatX_};)oXU|GgVeB+c66s~Vm;aJuP-!%>{iNhEN@$GDD=pR2|@=EZj`=ZBDo zP#E9Gtvho8xjLoe{alHVRigkVp=v7xtcWu|t5!>g^NU3GOFqVC+|DoKL&n3J+AeXV z&eJS?l#dVKftke*T@=s>pR&zW4ZozZao zlbWE=;pk^yR`lqvyQ4oe>K5=-x1ajLZ@-NB%oyrd2*U_4=`Tw#5u0F&B2kBK(l`X@ z0afFPb$AHlD|BLxL~;M;YgW!XTKGejtu`NIGN12kEGjWm}D=(5dmmv z1qUv;g^e9$GUN^t!@4MY8i;tKJaJl}lii6zh-VuG8Lf&S7I_soaST!8Ohq(K=XGNI zt>w^e$A_s9+8!URST<|LD9ExaV&s0ubO@v7DjN#RZLRnO)u&g)#~qP!2;U~Zg`tg+ zwc-_&|7=CP+&j4u!mEcU+euq=tq28|>~t;Ll?Q>4GtQfJ0b{YtvQ~WlrrPpDhZ``n zA$-Az4fbKyik91W&Ohg4*=BO5=2i%8z*nvi){0PY{fiY5z7JpLLI{H?BHqNcViaWm zaz%`f(B<*>x%i14zpls?wS{9+ zMDP4|6AiLNXh}mq%H3&6gVcFXG<<mE|e$Mge<4$KJc9gF6JA-vQv37#YL#-7PV4F}(yM+|OFIGq#`# zwJsQOqo5nhAC;&=#Lm;a2t_Z@u)$J+mcjChpD@2KM2jAT^mw$^z-uGsb>bE+Nn)J; z$;Sd87ySqAR*5S7{ty+q4@d>H^1g4LBQ)tWeo>h#6T#TRw?jl|RZSvRQq`+wh|m&f z?v-iT0a|__L`%Ssmg;-eM6jR}X!DhcV4TZGLqzDQNrtpi>1(4Q(C{nMFbo<#7NS8b zBhs*fnqM0!f!1G{l+7UJ;~`RlK541?S51UW+rSc7nTR7G;)g;+1oO4EIj}Y=0_$L9 zDjo$Dp9oPA>Xqe9gtd_pSPUzZ!cgKL4w1563~Uh{T=2aI;UWDJm=Plm+I~JM?z!-!=trL!8QsV z{YWS*-`=E#iBQAf*TS!>Wf<=cilP^rgfTGTdREt=Ya<~J-K})!D&qqHg}>(02z4B{ zw1XG=XakDfmlm1g_Cyvvy(1A=9EK`^2NJWWmKaA16>aF4dDG8YU5tqcz0reBYH_#4%F zTAvsI522~gqI?rHV)YsM7+!>t56zpxLa-9I{gA1`Oz+NM=XX&M)&_lL!b<1&m|CP; z$-{rjmiatPE#jf}-Urc;i@6xnA)KUfHLSEz?#GFc6~M{|TB2)eeiJjPc2>ts1hwM$ zqU&+^rr5eIG)9}*YKn^<+|`2<+XAf^`%V*M^>Ja!9@wre`T3J*TPnVopE4?GuIh6D z3rsEQ#ofsJWQeJOE^cKeFPhfgzy!CNrnNsLqst6QitMZSm1tWTlD*!LkbNq)f$B{fi&q1* z`M6N8E|Uv4j>Su-aiXaUr)PccL(QuagMbI z(d*HAd|Xao)?@S5jbxpg1u^JcQ8mnwFH;&A&OhLTdP45ya5!_E*ZJV!7A;cyzV4?kYTlC-xw48N;yonAIb@dbX5^o3%b)9?0G@H-?rd z3+ghNw=~KFWr4A0hJ@O?W@l4P)}qd?@RutP+}Z!qYZ2brKhOL5F&8jjc+G3RAl156-vt!;<(>p%~BqK2icF>0l zK7)*l&;Z=@(mT!xFWf<7TqA-7o@jIkJ*pmuB1(_`PNtgmLZkC2kHg=dagP6rCseUH z6Sk8=8}uRB>eYg{ev4q)T`dixefaRpA%PmKibH^#xXkIx!6y9^*tfV#phsb*qmN)G zZQ;fQ>{}Mrr~FOv0<@&zK@|EAJ`LWEj@bAJxA^9eJAVwI*nt@!0SkOIb5@Q}elCP` z(ZXT{wA|H@R+`dgOWG`&MYdzyRyWKizZ7j*I(zSdN=}Ui~nyn z{VX4Ug2%h&F$h*JqqJE~`SP9)+%6EnBF@r*`(Me~SCBBcGy4t141B|~qF=Wx>+4n| z@z++w`fDqi_?i{9zGlVx{>qA3e`Uoty^8>Ycdd@ZUs@g3Us{R8Uswt2FRae!pIe>Q zpIgc3S8*rLSFNt-pIP6i@m?pz|FO@Ko&sv!Q`}qye}sqdUN-QGDRCHoT8;>MmH!|# zBDAFTJ*-p{OK|z4qSwa-j3(_4VEj#Dw8UMlkK-;NV==pE9)EFJv~3*K>`FF~NW0L` z?YLD{x*oftBnjBX{60d&3w$iA_E9lHuY)S z{!s|wh9-nd%-AMz%Q;JobV=iY;%257;7ow&g{ij?Wr|1s^kzZKsiV9z&&Rj%G0n#< zJ{I_B=0pQ~(kfrfWd09{qZNkRqYyq&(`ym;5`LJ za)3(CJ?S&>T#w5MAvt93^z)T=h#aw&GsIf(KEym3?mL*`%2N+Rnf5Xy9%Y^=2{}$L zl3^$FY(BUw)2Hi2f$rlgwUs`APw}nK<8#LAQ7(fmSi^jo9|As0>Uri}#)Cak`2{?1 z+18HahYkLFJ`sVuh*+?l-s3MKmEIE!*58U>T=caBPVb^TOWKn|Q<{D_cAv&?g)=4i zU-&>UTh8Z^b@82u^A-Gc{)~^W^3jS{Ec&M$l=y}5BURUgtAu(NoxfZTK0kjs1fTDz zEzX~_d;t$c0U{jf4qv~WkGp)_=7R+B>O&q`)@FEMmH#vg{2U*@ipNL<|NjbY0>&#i z_MBG;D@5#HMPK39dFFO5KO^@Vv9&r#ss%3^iKBY59Vc%qAH^fO(F>d7-);o$;SwVT d^o2v!f}pi;b?S5R!_ggsR*xROeB!pb{|DywX$$}W literal 0 HcmV?d00001 diff --git a/dependences/__pycache__/pybam.cpython-39.pyc b/dependences/__pycache__/pybam.cpython-39.pyc new file mode 100644 index 0000000000000000000000000000000000000000..ca122ff88a06454c451666e45f94b6f4e9b03018 GIT binary patch literal 39119 zcmd6Q3wT>ccHX@g4}u^?Q4;mEEnm^LNLUm^JuJ&yNozeVdM#O&WP44Siv=+t1qlSG z3s4eaa5s&dSZT7Un^0N!R;r(r40S(=T81Xqu+kcGITY21(O=O}bx{rg?3i zj+^cOpSdrPAVID-`AV6WbLY;SbLPyQ|4c$JG;Rw;E7`Dy5!Od(@EHtI}$p8dm$&0d-IvQctPF z>S=XEJ)=g{htyGZOdVGrR?n&v>ZCfQPOCHO>^ld?den1c>+|c?IW@Y_tDaXMdDF;m z$oDL4yl>n!PviH#xz$kTwi-_xW9y$Y)C=nTn}#}n*BaY|w+re;yuFz3yBitXjQ5M` zCA_~R?^~30*+^Z={)J^2$%}In*=c*yDNozE^7M3}YENVj9Uh+CkJr+zyi?69d$R1< z7w4+i%B6$$#GHLDTPo!%_6xHGl`ZCHX41)1vyPK5Rf}_WxrB$BD`lq(IeRASR0^dj zyHK$wXG^(ipdeT#eN%<%wb_X@uv0RxcR0_qy4f}=cgq^KQ?)mH;`(iO$ z%H6O>OOxgEKRQseCD4E#xXg_E{$nCOT8;aioVl0+ktK8j%7!U8fAP{~uVgwoDovQP zPpHC7xsta}y5z9tse+TsszSD8U%G~oH*FV6HilFlaElm-vs2e>4An|LKV2af)K=NK ziVcu6c~qY+SA;Rsptpb_<8lLIZ3YZ*?2)~xjpqG6tBn7M{|s$o_d=H-T2Xqcl7bMgU1l<@%5%6JgP zqt<)u#HsVoUAp|rl~-Ro_mR^lFVs5EkDk78`P{{0AFVxMUjch^dg`EAg{jgs`x7m_ zeBwMO-E_8^uKk96qF5A8PRv!o@#-}+W;%N#Z_mz1_EHFg{Z2ls`cp$mcX%%I2jc;t zfo}q{&6%^W&w^_hd*CdGAU?1IPm-8IF`u?yuH+|Yi*113zJ^ho0|8Y|{j6Qe-}XUp zLPIgMTGkQAZjd&kGy0%}+ zmOZ(lIW>_74>2(Em{g<%BB7F<&fJ7>%BO3G?W+ecQd*!q`@)Im?9%M?1O|#bxk9Lh z*MMrWn4L;za}mSPeSHmLwk$Gi^t3&Q*)I$RQ5c~zj&0JmE0u1;<5K}3R;bunw6ju0 zYgM?{u3;3Q4TJ2xOwO-UX+JO^Jj)cPQ~A43daiP~rr9@wYP{bbgvOc0q?SD5O0@xf5>lFKiZeXx6e{Hs1lg=;AC6&x zOFc9_o4aODW-Al|GLxiC8{mYH%~fZikpeDCSAn-cn!l5shDI8KO2XJ@V`R92S(H=R zHqa|&QL^4N6Tw+5mv0aTp#vy-w+Nawd+c4uKzgB~-`X+eLKy;IS}Meq@>7_YP@2-F zK;~3zjz$^gdF2W3*$T8`z76@)6dJ5NlXpPww7tqsZUgoAyZ8Ay5XS(*k?5l^W*ILP>{d*^}70hnkX4ooA1e494Zqo+gVUvE*+^pEe zd_`1C$=)x*BryIk?fZxOITG6d?=O|1gsa(p8Nr;nAa=ifo6O1CP70X&{eBS}eMp8?s>+QaEy+p||J4G84krC8eF@eYpe!OUos;ZV2Mf zolr6?HOd?L8PTJQzO_@w=mO~f0vb4HPZe%){!_qm%!>{~_v^7AsAw(o)3voICD8q3 zx^~jO`c&(oc0M~ZL)9u{i#LQS5*6Hlb&2;Sz`|9D=SWk^Ucj6Ls7!2tt z$REsZkHMfY}m;_)K1f83l2Z)dLt7Scna68VH8;9}~3#OB^rc4A|)#iNac1 zMML&AjCdGndg&>$whe5(K-ZF(ii@yeX#oP4uo@Ru)*zRQ5Xd8}3uKn`11l!c;WZ34 zc89i7(zV~SuO4aLA8l3JDm#Sw56!{q=&RCjZTLnb!Pjtlva=iXkZYyn##XK5y!M4= z*#CFBPG;7cv|Ss1(x@G>Td2Ib@~jNX)~e3c;OW|dHH<^)q~Ub!ar^2st;gV%7tUPJ z5^8N)ty;{N(zS2J>UFp^;LB%TawXs3=?f>|u(&y!FVWkiXPz%^+b;Q9DYp}+p1X1$ z-UZjWAR6cx=IFjP@CFp6YjbG8{?-k6>E#nIcrCc;S_G9GR+?Hvvr*JoZF-yI`OHhL zR_<%DT!=rN-N5B^YfXMoQ|cNIqRAnEtGCO0i+Om{9Qz7Z_M@k%k=jvQotuGkhTb?x zM^WC|842Y@kx4_>TWS-i9<77aM3kHGLgm^i5K0u4uVxEHEXaN9K?8D!2Ox*l2@%UO zQCvf9Sbhy!+HV95e44HOIWgqH6W4k{zYwUCEluU^z55S5O)F&&{Pd(L+6QfOP(p70-2^J2ZGw&2&0!d%MsDp~lV!2k<&MNeVHi6qPLwf@H3;kT9 z7vKVbR$Yb(2lgL&`p}VshYx991h4H{L3L2gPE|5uW6H7zyv53nPAGS^NM;9v%H_Vd z@?QKuwK@iq@6!XVdBT8+M_c=LdA11W zF0*FQv2B2>+6vc;#XMvf7x7$lz)grP1Dftnr4M?N+MCtFc??m^ zn?j`xqx4+HQ1eZAwKn05wx+XtJVPS1;&45|t=z{|xP~kA4F{vE4ViSu!MR?7yNBZw zW59uDuo~#n5Qr`qS%f>T}moPxGGEl(_Y8XR)<_ zp?#xNzTJjWj_zr!`9!hnb+kK{XwA@#0#6PF5xdKS3ArkpEUR4LiuYwu%VGyc5w`HK z;B0c&=O;0ftB@6Iaqbs0V$l5ZA{H%Lll$D#`ICiw5$m>>WqX5MYwz_WJJI=?pIn7$ zjxE6>8i0+}{wGyNGr_txJDq5D-}PG56NC-}@F`!lM;R7`wudXCZ`ro21D(%3#oV>o(hYdW#0Lbf!=(tv7+jpI zXk2!hgOnXB;FprPZ>9h%9h=g6=iPN`xx1O}?w32jJtGUm<*=KNak$72)J~bx;tf^N zdZI$EqgD7Xi?9J5;(%5IX;5}A#de$K-pFJ3x)}x5ETW3UihzZe4lQ~D;FxXR=xSH9Br$}RujyiA`GB*4Bg!-gh!VCqDkvWkF*?0yp4Z_^KhXoZ@9Ck>T3?B-3k0W7N_i4zbm z70~|$#6w^pf*0hw5DG)BB!MBEt<&uc{u5AX1X}t$=DeQ@@d0L0aKn}SWV!6bUle0k<12}niX!hLOh_4(v>Q`dXYmbD(0{YNG9cz^0(BZ4+GBG`g&bkYi!NQg!!-rZz^xeQup zB$z2FDT$#ydRkY))Sk}kN;+36aWkyprG^_Ftd2aL#l0U5QBu9c65hgS^aMFvJZtDu zLVpSz$sP^xEH`?HB=jjsDB2<02-3qEn&fM8=(N-1NUDR?U^Qx@nUHqS*`!NJY8^}Y z%lpQ9|5bON3lGcxhQNcA(eRj%<^@X)O?c>CfQ3x_#(K6aoFSbcwykcBP72(>#&a3s z+kA~a4#rGi(8(gAu~BkG49;cr6lkp1d0h1ALR%#xP%{}<-!*y{eW)}xNKPQ?%1)z4 zXGvThq+6p)-d$Xs^1@M+{dW;WAfC^{UdSQ7p)%s?aVq++DunK;ecBB@yok_42<9_7 zGQ%Fn+@fE==LIHU2-)DkzUjGrG&0iQ{6Nav|I=Z8ekBP}!2n4nuiBTZ*xkz^q_8C( z6dK}KUlGmF)d6?t4Gc+&H3r6g?0QRvK4VWZXnI(TB%CkYfedd!fsVkik1gNAEQ3lw62yX>q)5ZL#awu2EHKpk1Og?eA?%OC;k&Xyf)FpRoZm@PIQ zA~r*Y8fQaWel2E|XTc2!8)Vekua+IeehpsBqBF)DX0l=AVA}XUnB7JPeox`I@+=-x zM%@@c^13x|&Rca;8M8*^Vf-~`4Yb^8QXxZYvVqpS-zxYx*n0{s-0-5rA8GD+Y%2)PH zVwJd0IF%skUp<2Ncu^nUsM+*y$9)+)Fnv|A$?h0 zk@Q}qUskV3I*s(J>Z6k0hxBV|Owz+hUsbP5dOySK~VfOJNET+#=T&Z-GXA5uA` zFfN`_dCb76MwIF~JAwdbkvJ>r+^Blk&{s=K z)vu=;MgMA2mes+llCaAd0K(TMt=L)*FKz82O}u7REP>)pm_;A!wys$xuYw0(6Ro@| z@>*%;)&%~Q>1Ivm>hyDK(wgnmMaE26LlafM?npfYa4H2J6TP%{SCa151sspu@zVPm zShH0M6JlXGl`Iq5{iOFMVL{)*8I!8Lx-%_ea&y65#Qq2zpM}z7SvytE=6F8B``Ls9 z?F8zq%9rp|<~=I|YXMJK2Y6ZqrXzJrS+l0|^VMkGx@($7HO9Zy_^erZw%Wn8*b-ya zgo;!<_1C*RCIy?~uKD$5X{K(SHpZ{dN9)mrq>9})E4L&Szi&F9u1Bl!YS%(H4_MtY zou31gj%sH;Di9J|4Wxkq&^k4=n;%=YMA9p9&v4$a#_j=PwG&p(AKf!l7h3SX>b_?H z$NPY=PAE$t)x*?V##`oF)?`FhWsNJ?lE^ zKW)s%-!^YFN~JEeqRt|tWZj{tM0hug>gZ2<((yP*kwxb`2lto?q#DW z8~^qBu4N!v&dEibux|DP+SY&`w?*A9kt9s}@Ig_66srO(sKYZVOJK_x0*MaIDe0}7;skUJR?9{bM_1lq~_)feS z6%psn>NfQ#p}|k!G{qiE%&)Jn=NLrpc5Ri-_tq`I`a-?; zo`v~g-NPK2*EMLVCFVDlx~2YreyeXp-5zDYmKeO?aJg`{b}^Ye%RM+avms(AG4#4q zP7q=lFJ}zkal<~-tpaz5;Zj8SM!B>%KRtsGWg5k4`!(!!PVqos))qTmmI)4pOKod1 z=`PsyN*1d~ptJjHJqV@Ht4cv%yG~a=_cpmtDvK=HxIlPpewjau zpLz^~q%u3<9&~hfIkMt2pvxfGDoSQQi=mfHz6`g<-pZ^Sq8(m~`-9G-5#tj0^%0+0 z3iVk1odmMb~3fjSdSQXc9jxLcKUp#K$m_%%`WYW3qx5g3W?IMTQx9 z8O#9jfp8)O9F+q|2sV0So<+la&_e{|lMCr2h>@Ux6)-Q0R+^)YDWJyI3-%2$Ei{t_ zJ|ef=4MBHn3V6e{JkK#Z91aA2HrDHD6uC_Qc({8x>`LWn}f#jbaREri#M zk~M79wp7gp$m8p=h0go%McjwS;#K43@B4M*O?Ba)pyqxJC5DZXa`7jD|5xu(QpP?I zygtyH06w%SuszUPSSs?QJVz_2mnc~!Z+(j|z3{ChZ&r~KJ!1Sf+=IdMT8PHw5_ssj2+JN_Mbwh3Z6nZMNy|kO()uuZ6&R2wnkqYrU;#?u^C522b;NFkAUkBgX896eph!c zJW6Xr#>n;s8#Bg0j}vcGAJhpK{jzwhiT&_9tL%x$&3^@`F+IM*<4mn5^g8AMz9!U8 z8QD>cY1r@hvH`YJOg%auS2k>98#eE*>;15CqZpa*=u-8LZLrQ^;YOXWuSe@0+HPuw z=Ti@;Ebq=1z+?^7iEY)knpp?uEAPA=zxj^LfC!*|4l{XI^|891$^BYt)GpRLs*kJ3 z1molIcy_2KGz?g+_I8YM^$6x@vL0J_q7HiX_j*$8);6gcxNC_;nw;-=73+f1q#g}* zQ_6xK)ELwd0jEn1&UXvGE_lp@R*h$VUACc?QF|~}V$NxS)$7u-;qBPX(>@$P z&@>QM9oA}n9eQd=x1Vdd)?Zuas@*O<8<+o%-bCwxXCpYa$>sbewO4c6on3MD$Pui& zV6WRU+K;I8!tOfOl$4gV4FAe{0A5=KAJ^R9)63uVU^+@$YWLG)f!y8*g`1b`FCc{(B9- zq50q8e3UbMvAm@e@FlL60mF4j+ocj$O0HGrb~s;Lek_DnK4 zDq$I7!0=WC0`yHn7*ej`R5qG&31H1>ld~ko~ z0jxTs@CRS{%2ytA^Yt&8dLNH5>*T2iT-80;hlleL9*x*^cB+u8b-Y3T*!`&|8&P(8 z!$Jr`hsxh7&4RU29)~~-MAd1;GVITWzWwHm<6zG^N^9M@L}+n-ER`EkDHNs6%>8|w z&P9<%0)wQ;Jm)kEcOWoDs-7e@$C(vDZRZ42v0M=Wv4}>Pop56$8?lQ*S%;IKDQ0u| zM)cz7#WRhrOi)C+qLI``24rt3dIKRRjgA*DyQK0?6ZC#DpPy+&XUc3_4I4_0WMDhf zO+Vu}5mFK7sfm(mM6Z_-jgUaJ)Qq&OWBP{TOC!eHX7Z^mjgHf2UO|NN3(gelMY$v1 zh+zGOMrUtl_&V0(E&ejOEDq<7ZIVC+oO;MlRx{Od2G%;J;6!CCDM$BbfU@em32rO_>fK?v!*BCjWoIgRv}iZN1Dr6b=kOcUgUcxuf)M2-m>#Uc&%XEW9)ZF zy3LK|NWe!A#)J5dyyYsTqhDoe$tG(|A=iPU~IuV^LGN$Zqr7Jltq!>gTJJPc)HDp z0r?1?-DYnz-kCtTZj|C*olaPN(T(O|a~S1XP;AvCz%!l66iS6^CX+rYo0*m4X`fXx zA4Re{uYzB-bxi`ahAnI_1~@p%gGGDRk9$ieoj5C6M~fnGa6L~YS13Hd%>PQwdc9Oj zJbU@2*Q=K=jNq;MfCD`H4!Y(nSJ=EB3eywUf5=i&q93&gk^=?$fI|?~A}e&Rpm|2j zW?RHJ;~&J~a#Tl_RWUGaJ#6jUAQur}^bm<#qgcL$faF_+{B3c{Fyev7b#<^9v+cd- z^Z~Cdln3JlB?1p;nsgA33s-i1Y&$rhBz_Yvig7%vLKnKUDMY=*U)h7+L8cTMg5vPQ z_*cEym*n~!aa!;w00Q^nRkHUUCxK}Z-4%dE<8|BMHNdSK@Qh`Pw|SoudWJkeVK;ys zorYj`x#|la^d$;t)QBhn%_^5myuw5uI~%cUN8lo&8-&1`8jZt3-l-Fw-j~R9pLzrY zEr}uu$ecemDT;_M*JPMzT0dVD)GHy6RS?ib(pJ? zQ#!9p^1c<^=WYP_j&JMSZ$0wuGx)aN|JLh%+rV#?@4=0Bm#CcjLdTV^v_0=pXJ zen)L}|8C=7_&ECN5$6v}39lahZV&3K9aFNFQ=7quhjfj~pnBM?JAu`z=?s^))w!qh ztXEf{JwhnR+37(;`yX}NzZo30U1+z|E*ILruEzl76Qx0cN>22%t_%Hp@cx**lh(&w z=ufnS9)~|kT37icwOdk_(_hDql=G{lBZ391dnI*UZ13ZJz{h4e-FP4Lv87sIdD4Z; zR`2n?t09-;d)==e!`HO`b)WuPFtFku#%garMhI7u2aYllsY4xf{~kJOyk%gAgOMKz z?4BiXX`@`wSns$S<6aM9VE8`bzxU~j8~yj=^8R7}-M!4wf8U{ax#Ma7{fxYy_1~Wp z_~-oheeyo)zxT@f^Zxrs3Is*Cj{hIUxV2p6N*a6MVbDZcMoQ`V9&?^_=X(rus8ikaOF2?X&*PYL z(VfRJ=Nw8^{8CjZrRQS|{>@G`D>)(qK|>WewVUy^q%%#@`1#CGoNee!~H7Q+3k z-hE5RAy60|qMT}*WGaU9G2)9oiS`5#cB0?5yp^PWZ2$g?@u`o6CCl> zKHLmMF2raC!>q-Cbe@13ASh|yX2RR+!-JAlnewEL!3W5Lhk=`9gd)J5!>Dsbb72<& z^%KU2mcWppdBA9fxyK$n^dun$!Qv@rgkXemLxgp5Sy-*;0fkFJwx9>0SRraOb*Q>` zaB#0~Aqu4S?K_Z4!DDfR6bV^ut?th>X*;%rHo~=63AuKuZFM;b(Xvgn1ZihKd33H_t(BM`23KQvu>SE!Q8WHEbE^q_a8P;?_z%js zFcuLXwHN)CaIlG+25cO7a=@GVdIc5CeI8b)Yb1PI=2fZ712E$bjbM>Iq5jvPxjt+H z&6%$+KMKi~W#I><$DPO6od|Qlwy%WU<5rR4q3zH&8RUce>>>xmk|HmgcV4(dt{0Fa@pYU<4)8Anh3Y9Pmfltgrl&84}W!7`S`Chcp`G8>n%QteB8n#6?4k?(})S*8*#Y1 za2-J-Lib&yI*nx-Jy&jII8HH!)@vP%AJ5Fr{8q49qsylo>$6@wui)7WH&4Q>9j2SR zg%gKyA^bYP*XpkFU_HYxM7}i=+1RAe0-b_ zYW4@B11akta5_DpihW??A*%d=Xm&KgaGt`W_MFdU*bgCtL8QOW-C<8+$=?zCP}4hJ z#>IA61WnI=x9BlFF@WW2bL|oP$WmyJETs*0MNlpTt-1Dyz5hT^Oj6^p?3T1Nr2DDo z{FXM?2DF1qpgp=2njE_E3WlLI)dsZYrju_6mqLRF$Ai^OTT^X7I~){~)VTAZrO*PK z1ZBdsHPr^RW=4{>txKWt%%YDq46Uj5h<%`W$Pn5yOQFd*I}cn4T61kcYicX?0X(8jxc{5WOq5YxDsl$%zStIXeR4^QmJ&v;5gJDJ6e zLE0>X1ET}D6~zOH6J)hFn;NHK$!%I69FN~zud_7kg60Y~kU3tPZ$X6{jR*^GAkAI0 z!K+@c;L9MJMGmB3f@-Q;aKL8(Yc9ux@l@@dzjpMi=S`~ii6-he;sc)GXpSmn-QGt^ zcrp&2H0Szb0(@X#B1<;&Smj` z<1X7KyfFj*Ed9Y|4T3|}APeKVWe%Z(2UA`-^2hzcF%$Zcxz2SNhRk>U^`TwD7IuOD z3FeZPU*qGc>?SM}nuzm{&|}O)Q=4WXT>xJY?nqhB1Ax1 zuA>_qczt+a$OU~cwP(NIJg+jLb4NQ=^{Q9irJj%>+MBB8i($$zN={YF7w6j0=E2=F zGk0ia?yA?YX2>3YT%z3Njerko1Zx}*W{6hI-r;7L&D7%Rw}F%G_~^ntvS-kxVBlzB zUacK|e16czXU*pE~ znM;>myflK?7=`F|y_?J{udstF>o3_5N4MzhXJGjF`=F(J<-% zc-ItF;2FUE^4?u`qqxM4_eWQ<`bI+=cdbY~H=x^UmbbAZAp2InSf1hkekDDygm-?TDW zxP$6EM=&Gmg7nVu0KiCUIcl2uquI5%_YRI=Dzz93xHqyo*vwO^ANGXUJZT+FGbsG0 zi8SPOgd{jsQq9%%$it-82~OeKPxyTFfd?i=r)y)JfG!vq8NhRJpddc0lT-o&xl@(P zqXXJqeG;m|FEw~7O;JYUgPB9Dil??%um&6ij-gZVwbK8^l7pvG$Bt=tl6HQ>n8A19 z^NFF8L#K#OIOmf|dn}dZ|H2{M1d1$hNgGu_+DX(jjhM?d8TFc*iyZBW^K;yC^vs9K zxw%}s2DyVF@0~WaJigHm30(N$Z{gaC~QCqaP@^d_doW5SQ>Whm(nuGkIeKXx}eRNEy?S_ zTOT5@Ec_;y<}9DIdHG(heJ$U%}FG1qQ|(7XrIf*yXJ*=E$6$Dhd3DDqpdq= z0l7M*gZbPUo>k-i_+KuV%`duII-H*-vY+uWHsf}F79TPm*3=FO5_O(s>7#sn7!S-A z9w>BfFde0g6?d>xVgY)G5kukH=ViCd0i4L>?mmbgUP05I$>8I2~#<_K3b+lsOQ(xDp$enou~(RrvTH=w-&@ zMKRjPwg#bV6@lqpPz9DC7V|z){uEKJ12=+;iG;(tv;Pi@NRKDHC`A%l`2-#YW{3Qf z8MV#afxms24gFTX6~zqPVJ0HIpEMp#V8$?7Gr_rlc+V~P_Be9*++xwC&cG)=cgl0O z#VE{f)asLW#92lVC)$IV`$^*^#!>d5R&8C=*l!Z<6c-%h?6Ha-<8^n8hsM?dzUo#| z%KNR9;hLF1{R&|i0Ve%r2_|9_Oi?82&?p**09`;<3mga!b~2ef;A_?4w?*=o#}F3` z*F795Qy+@R}+S6DoUP1YfSH#P0kDDR9dWf=}v>?}t zP;kk{)v}B|2!xz*-mD84iw%;s;uBOKU$KGQJDCmP3r1|Pg|b$(g7KfOh?biyw?k+H zzS6E=D?-8b&sIeE48G2V5C&63e1L1kD9HZ#iWr~3*Ed2KL!9+@@z;u0F#d}b(Q;E~ zK7=-CsH`#8icygL@QN6@@$<CEEl7@biyW5fmsq=wo_*K%- z+k8(0D>6LADYe0PY5V7|6C2i8VKU>&SX#bco2oe&kFURmBmSQ{yU#jr9d3?%-Z z5Gm`$z!uTL1>Xk{9?~y?8L=`M47^?lkr51zmbL_ifh^Z3e<9C(sqS8_W!=ZZjf*I; zVAi8hShG8yM@CCQ+J_I{8zLhZ9<6N(wo&NlM?zuw_9it{tXA$6~#3 z5g)!ENc1v$!0%cyPV@$Zf9Nvgk7f8BA__@eU=BI7M_Mv%7e0I+3VTbJU^EAKnwrg_ zz<5vL9PL(RE*J_L8QjB_y&mWAH>!2CK0E;aK}(-S`6g&hOne+)g^3T(o5DV@5!d&S zp~6V-@?YoIQ4p2}y01pOZIt@~ z6!zG-xh1-$rZ+K?>SlF}M1U&JCAyx5uZpGHLu0d9(SAi=4osq7EM=rTi* zVp`@iMUxE0UT-M!ju!XsY4ZVR4K9C=U`AU7>uu8GTs}v~A!4)SlSF#gQ zy!0hruZalpAJH{%l_jir*1=JtIMd;}uALvpF}8T=v%uMbD^^TrC+a2;*CDz<~_K8?kz0or_AfK``Cg=@p&rPDaj)P=)LxY~z@nef(q*Ob4wVv~RM zH#~j#)o#J7%dfkOqq{hw+Jgx6XgxlOgQ|Fs&09BLn5dchljiLbXfBy^re2upIai%`J*FTa5fk|OF|CWz$;N#PHaBr|smEDh2hr2s)T@FP`K30RIO$5xJl@u!F$68Mes2!4@2lPH<-dH!F;ok6S&m z|7J^xxLMoi%LCZ~=625#CE?pkG{OTVfwAWv%oE$>J+rf^CTmeoSNQueJaFOuE3Zd* z;r=`?*T-DIeBpuWoxE^guGdE%^lpSY7c={XdIT5gMR3+4vVi%IOLQZ=_Gqh7ic1

iYM|r~i&Wv;X=X`-0*hQf&^pITjYCufCMUd>OmWI$UKKyJ*mbbE`kSVpQTc-NO+GYSgB*J6^4?IR+poXQ zPxWvF()Rw}Z2DO~{s52n%wrI$Trz1hnzH2sowQvbfJL0$h`j$y&VCdLgL|!CK@`AO zEGzn#mSz2=6-oSGD`Nd$E1G!Uidyemv5kLW#jL-u;(hNSWZ*rkBk|`}hxO-HBJpQd z!um6-Gy11ir}d{+GWsXDPUlaouIL|IU#amLCdL0F&yt=3YTZ-ZNCkg{hwxN3@QNvM z7=K!h2zrozGc+Q!l=i);R1r&X`JU(z0mV zIHlQ@Y$B0%p`kl)Evj@qc0frIu#4G!goqdTSkzgxHGeCFy;s=kHM*76?X$NVRXyC& z+@D4U5Uyj2$z_!P2!ewmKf=N#sR}EOfSH305c0yZz0MQkL>AI&zU2guvze)WK4^hBX+HMzA!1kbjSTjk%;Q>94?>yrG9Vsf zo~Q^pMK5w;7xQdBxERxiutj|i^Od?vpS$1A6gx!k=yUSup(fqWUUt6BN5ExCy~w;v zc(5NTKZOTw$=aFxq``kLCL)j)5ert+2mD3k(FcOT`rGh}JGz#@>0Oj(Nqcf=O4AQV z?z8x(}vdkB>WipbGAlgF6L>jkWO|*y2CR0zb~j&*L!` z!T)+f8-DQ$4lL&t!u}AOSJ8Lvb)LDC%g@LKMQo`Kcxu6iM&hWRXvfLd%D3YY-Q?{aC`&y?yA3q%3IcW9hq06VyYyLm|k8Gj< literal 0 HcmV?d00001 diff --git a/dependences/pybam.py b/dependences/pybam.py new file mode 100644 index 0000000..9f02c02 --- /dev/null +++ b/dependences/pybam.py @@ -0,0 +1,755 @@ +''' +Pybam from commit ba460f1 converted for Python3 by Hannes Luidalepp. +Currently only dynamic parsing is functional. +http://github.com/luidale/pybam + +Original description: + +Awesome people who have directly contributed to the project: +Jon Palmer - Bug finder & advice on project direction +Mahmut Uludag - Bug finder + +Help: print pybam.wat +Github: http://github.com/JohnLonginotto/pybam + +This code was written by John Longinotto, a PhD student of the Pospisilik Lab at the Max Planck Institute of Immunbiology & Epigenetics, Freiburg. +My PhD is funded by the Deutsches Epigenom Programm (DEEP), and the Max Planck IMPRS Program. +I study Adipose Biology and Circadian Rhythm in mice, although it seems these days I spend most of my time at the computer :-) +''' + +import os +import sys +import zlib +import time +import tempfile +import subprocess +from array import array +from struct import unpack + +CtoPy = { 'A':'"' +} + +wat = ''' +Main class: pybam.read +Github: http://github.com/JohnLonginotto/pybam + +[ Dynamic Parser Example ] + for alignment in pybam.read('/my/data.bam'): + print alignment.sam_seq + +[ Static Parser Example ] + for seq,mapq in pybam.read('/my/data.bam',['sam_seq','sam_mapq']): + print seq + print mapq + +[ Mixed Parser Example ] + my_bam = pybam.read('/my/data.bam',['sam_seq','sam_mapq']) + print my_bam._static_parser_code + for seq,mapq in my_bam: + if seq.startswith('ACGT') and mapq > 10: + print my_bam.sam + +[ Custom Decompressor (from file path) Example ] + my_bam = pybam.read('/my/data.bam.lzma',decompressor='lzma --decompress --stdout /my/data.bam.lzma') + +[ Custom Decompressor (from file object) Example ] + my_bam = pybam.read(sys.stdin,decompressor='lzma --decompress --stdout') # data given to lzma via stdin + +[ Force Internal bgzip Decompressor ] + my_bam = pybam.read('/my/data.bam',decompressor='internal') + +[ Parse Words (hah) ]''' +wat += '\n'+''.join([('\n===============================================================================================\n\n ' if code is 'file_alignments_read' or code is 'sam' else ' ')+(code+' ').ljust(25,'-')+description+'\n' for code,description in sorted(parse_codes.items())]) + '\n' + +class read(): + ''' + [ Dynamic Parser Example ] + for alignment in pybam.read('/my/data.bam'): + print alignment.sam_seq + + [ Static Parser Example ] + for seq,mapq in pybam.read('/my/data.bam',['sam_seq','sam_mapq']): + print seq + print mapq + + [ Mixed Parser Example ] + my_bam = pybam.read('/my/data.bam',['sam_seq','sam_mapq']) + print my_bam._static_parser_code + for seq,mapq in my_bam: + if seq.startswith('ACGT') and mapq > 10: + print my_bam.sam + + [ Custom Decompressor (from file path) Example ] + my_bam = pybam.read('/my/data.bam.lzma',decompressor='lzma --decompress --stdout /my/data.bam.lzma') + + [ Custom Decompressor (from file object) Example ] + my_bam = pybam.read(sys.stdin,decompressor='lzma --decompress --stdout') # data given to lzma via stdin + + [ Force Internal bgzip Decompressor ] + my_bam = pybam.read('/my/data.bam',decompressor='internal') + + "print pybam.wat" in the python terminal to see the possible parsable values, + or visit http://github.com/JohnLonginotto/pybam for the latest info. + ''' + + def __init__(self,f,fields=False,decompressor=False): + self.file_bytes_read = 0 + self.file_chromosomes = [] + self.file_alignments_read = 0 + self.file_chromosome_lengths = {} + + if fields is not False: + print(fields) + if type(fields) is not list or len(fields) is 0: + raise PybamError('\n\nFields for the static parser must be provided as a non-empty list. You gave a ' + str(type(fields)) + '\n') + else: + for field in fields: + if field.startswith('sam') or field.startswith('bam'): + if field not in list(parse_codes.keys()): + raise PybamError('\n\nStatic parser field "' + str(field) + '" from fields ' + str(fields) + ' is not known to this version of pybam!\nPrint "pybam.wat" to see available field names with explinations.\n') + else: + raise PybamError('\n\nStatic parser field "' + str(field) + '" from fields ' + str(fields) + ' does not start with "sam" or "bam" and thus is not an avaliable field for the static parsing.\nPrint "pybam.wat" in interactive python to see available field names with explinations.\n') + + if decompressor: + if type(decompressor) is str: + if decompressor is not 'internal' and '{}' not in decompressor: raise PybamError('\n\nWhen a custom decompressor is used and the input file is a string, the decompressor string must contain at least one occurence of "{}" to be substituted with a filepath by pybam.\n') + else: raise PybamError('\n\nUser-supplied decompressor must be a string that when run on the command line decompresses a named file (or stdin), to stdout:\ne.g. "lzma --decompress --stdout {}" if pybam is provided a path as input file, where {} is substituted for that path.\nor just "lzma --decompress --stdout" if pybam is provided a file object instead of a file path, as data from that file object will be piped via stdin to the decompression program.\n') + + ## First we make a generator that will return chunks of uncompressed data, regardless of how we choose to decompress: + def generator(): + DEVNULL = open(os.devnull, 'wb') + + # First we need to figure out what sort of file we have - whether it's gzip compressed, uncompressed, or something else entirely! + if type(f) is str: + try: self._file = open(f,'rb') + except: raise PybamError('\n\nCould not open "' + str(self._file.name) + '" for reading!\n') + try: magic = os.read(self._file.fileno(),4) + except: raise PybamError('\n\nCould not read from "' + str(self._file.name) + '"!\n') + elif type(f) is file: + self._file = f + try: magic = os.read(self._file.fileno(),4) + except: raise PybamError('\n\nCould not read from "' + str(self._file.name) + '"!\n') + else: raise PybamError('\n\nInput file was not a string or a file object. It was: "' + str(f) + '"\n') + + self.file_name = os.path.basename(os.path.realpath(self._file.name)) + self.file_directory = os.path.dirname(os.path.realpath(self._file.name)) + if magic == b'BAM\1': + # The user has passed us already unzipped BAM data! Job done :) + data = b'BAM\1' + self._file.read(35536) + self.file_bytes_read += len(data) + self.file_decompressor = 'None' + while data: + yield data + data = self._file.read(35536) + self.file_bytes_read += len(data) + self._file.close() + DEVNULL.close() + return + + elif magic == b"\x1f\x8b\x08\x04": # The user has passed us compressed gzip/bgzip data, which is typical for a BAM file + # use custom decompressor if provided: + if decompressor is not False and decompressor is not 'internal': + if type(f) is str: self._subprocess = subprocess.Popen( decompressor.replace('{}',f), shell=True, stdout=subprocess.PIPE, stderr=DEVNULL) + else: self._subprocess = subprocess.Popen('{ printf "'+magic+'"; cat; } | ' + decompressor, stdin=self._file, shell=True, stdout=subprocess.PIPE, stderr=DEVNULL) + self.file_decompressor = decompressor + data = self._subprocess.stdout.read(35536) + self.file_bytes_read += len(data) + while data: + yield data + data = self._subprocess.stdout.read(35536) + self.file_bytes_read += len(data) + self._file.close() + DEVNULL.close() + return + + # else look for pigz or gzip: + else: + try: + self._subprocess = subprocess.Popen(["pigz"],stdin=DEVNULL,stdout=DEVNULL,stderr=DEVNULL) + if self._subprocess.returncode is None: self._subprocess.kill() + use = 'pigz' + except OSError: + try: + self._subprocess = subprocess.Popen(["gzip"],stdin=DEVNULL,stdout=DEVNULL,stderr=DEVNULL) + if self._subprocess.returncode is None: self._subprocess.kill() + use = 'gzip' + except OSError: + use = 'internal' + if use is not 'internal' and decompressor is not 'internal': + if type(f) is str: self._subprocess = subprocess.Popen([ use , '--decompress','--stdout', f ], stdout=subprocess.PIPE, stderr=DEVNULL) + else: self._subprocess = subprocess.Popen('{ printf "'+magic+'"; cat; } | ' + use + ' --decompress --stdout', stdin=f, shell=True, stdout=subprocess.PIPE, stderr=DEVNULL) + time.sleep(1) + if self._subprocess.poll() == None: + data = self._subprocess.stdout.read(35536) + self.file_decompressor = use + self.file_bytes_read += len(data) + while data: + yield data + data = self._subprocess.stdout.read(35536) + self.file_bytes_read += len(data) + self._file.close() + DEVNULL.close() + return + + # Python's gzip module can't read from a stream that doesn't support seek(), and the zlib module cannot read the bgzip format without a lot of help: + self.file_decompressor = 'internal' + raw_data = magic + self._file.read(65536) + self.file_bytes_read = len(raw_data) + internal_cache = [] + blocks_left_to_grab = 50 + bs = 0 + checkpoint = 0 + decompress = zlib.decompress + while raw_data: + if len(raw_data) - bs < 35536: + raw_data = raw_data[bs:] + self._file.read(65536) + self.file_bytes_read += len(raw_data) - bs + bs = 0 + magic = raw_data[bs:bs+4] + if not magic: break # a child's heart + if magic != b"\x1f\x8b\x08\x04": raise PybamError('\n\nThe input file is not in a format I understand. First four bytes: ' + repr(magic) + '\n') + try: + more_bs = bs + unpack(" sam.\n\nThese two headers should always be the same, but apparently they are not:\nThe ASCII header looks like: ' + self.file_header + '\nWhile the binary header has the following chromosomes: ' + self.file_chromosomes + '\n') + + #HL - decoding header + self.file_header = self.file_header.decode() + + ## Variable parsing: + def new_entry(header_cache): + cache = header_cache # we keep a small cache of X bytes of decompressed BAM data, to smoothen out disk access. + p = 0 # where the next alignment/entry starts in the cache + while True: + try: + while len(cache) < p + 4: cache = cache[p:] + next(self._generator); p = 0 # Grab enough bytes to parse blocksize + self.sam_block_size = unpack('> 4 , cigar_codes[cig & 0b1111] ) for cig in array('I', self.bam[self._end_of_qname : self._end_of_cigar ])] + @property + def sam_cigar_string(self): return ''.join( [ str(cig >> 4) + cigar_codes[cig & 0b1111] for cig in array('I', self.bam[self._end_of_qname : self._end_of_cigar ])]) + @property + def sam_seq(self): return ''.join( [ dna_codes[dna >> 4] + dna_codes[dna & 0b1111] for dna in array('B', self.bam[self._end_of_cigar : self._end_of_seq ])])[:self.sam_l_seq] # As DNA is 4 bits packed 2-per-byte, there might be a trailing '0000', so we can either + @property + def sam_qual(self): + return ''.join( [ chr(quality + 33) for quality in self.bam[self._end_of_seq : self._end_of_qual ]]) + @property + def sam_tags_list(self): + result = [] + offset = self._end_of_qual + while offset != len(self.bam): + tag_name = self.bam[offset:offset+2].decode() + tag_type = chr(self.bam[offset+2]) + if tag_type == 'Z': + offset_end = self.bam.index(b'\x00',offset+3)+1 + tag_data = self.bam[offset+3:offset_end-1].decode() + elif tag_type in CtoPy: + offset_end = offset+3+py4py[tag_type] + tag_data = unpack(CtoPy[tag_type],self.bam[offset+3:offset_end])[0] + elif tag_type == 'B': + offset_end = offset+8+(unpack('"' +} + +wat = ''' +Main class: pybam.read +Github: http://github.com/JohnLonginotto/pybam + +[ Dynamic Parser Example ] + for alignment in pybam.read('/my/data.bam'): + print alignment.sam_seq + +[ Static Parser Example ] + for seq,mapq in pybam.read('/my/data.bam',['sam_seq','sam_mapq']): + print seq + print mapq + +[ Mixed Parser Example ] + my_bam = pybam.read('/my/data.bam',['sam_seq','sam_mapq']) + print my_bam._static_parser_code + for seq,mapq in my_bam: + if seq.startswith('ACGT') and mapq > 10: + print my_bam.sam + +[ Custom Decompressor (from file path) Example ] + my_bam = pybam.read('/my/data.bam.lzma',decompressor='lzma --decompress --stdout /my/data.bam.lzma') + +[ Custom Decompressor (from file object) Example ] + my_bam = pybam.read(sys.stdin,decompressor='lzma --decompress --stdout') # data given to lzma via stdin + +[ Force Internal bgzip Decompressor ] + my_bam = pybam.read('/my/data.bam',decompressor='internal') + +[ Parse Words (hah) ]''' +wat += '\n'+''.join([('\n===============================================================================================\n\n ' if code == 'file_alignments_read' or code == 'sam' else ' ')+(code+' ').ljust(25,'-')+description+'\n' for code,description in sorted(parse_codes.items())]) + '\n' + +class read: + ''' + [ Dynamic Parser Example ] + for alignment in pybam.read('/my/data.bam'): + print alignment.sam_seq + + [ Static Parser Example ] + for seq,mapq in pybam.read('/my/data.bam',['sam_seq','sam_mapq']): + print seq + print mapq + + [ Mixed Parser Example ] + my_bam = pybam.read('/my/data.bam',['sam_seq','sam_mapq']) + print my_bam._static_parser_code + for seq,mapq in my_bam: + if seq.startswith('ACGT') and mapq > 10: + print my_bam.sam + + [ Custom Decompressor (from file path) Example ] + my_bam = pybam.read('/my/data.bam.lzma',decompressor='lzma --decompress --stdout /my/data.bam.lzma') + + [ Custom Decompressor (from file object) Example ] + my_bam = pybam.read(sys.stdin,decompressor='lzma --decompress --stdout') # data given to lzma via stdin + + [ Force Internal bgzip Decompressor ] + my_bam = pybam.read('/my/data.bam',decompressor='internal') + + "print pybam.wat" in the python terminal to see the possible parsable values, + or v==it http://github.com/JohnLonginotto/pybam for the latest info. + ''' + + def __init__(self,f,fields=False,decompressor=False): + self.file_bytes_read = 0 + self.file_chromosomes = [] + self.file_alignments_read = 0 + self.file_chromosome_lengths = {} + + if fields != False: + if type(fields) != l==t or len(fields) == 0: + raise PybamError('\n\nFields for the static parser must be provided as a non-empty l==t. You gave a ' + str(type(fields)) + '\n') + else: + for field in fields: + if field.startswith('sam') or field.startswith('bam'): + if field not in parse_codes.keys(): + raise PybamError('\n\nStatic parser field "' + str(field) + '" from fields ' + str(fields) + ' != known to th== version of pybam!\nPrint "pybam.wat" to see available field names with explinations.\n') + else: + raise PybamError('\n\nStatic parser field "' + str(field) + '" from fields ' + str(fields) + ' does not start with "sam" or "bam" and thus != an avaliable field for the static parsing.\nPrint "pybam.wat" in interactive python to see available field names with explinations.\n') + + if decompressor: + if type(decompressor) == str: + if decompressor != 'internal' and '{}' not in decompressor: raise PybamError('\n\nWhen a custom decompressor == used and the input file == a string, the decompressor string must contain at least one occurence of "{}" to be substituted with a filepath by pybam.\n') + else: raise PybamError('\n\nUser-supplied decompressor must be a string that when run on the command line decompresses a named file (or stdin), to stdout:\ne.g. "lzma --decompress --stdout {}" if pybam == provided a path as input file, where {} == substituted for that path.\nor just "lzma --decompress --stdout" if pybam == provided a file object instead of a file path, as data from that file object will be piped via stdin to the decompression program.\n') + + ## First we make a generator that will return chunks of uncompressed data, regardless of how we choose to decompress: + def generator(): + DEVNULL = open(os.devnull, 'wb') + + # First we need to figure out what sort of file we have - whether it's gzip compressed, uncompressed, or something else entirely! + if type(f) == str: + try: self._file = open(f,'rb') + except: raise PybamError('\n\nCould not open "' + str(self._file.name) + '" for reading!\n') + try: magic = os.read(self._file.fileno(),4) + except: raise PybamError('\n\nCould not read from "' + str(self._file.name) + '"!\n') + elif type(f) == file: + self._file = f + try: magic = os.read(self._file.fileno(),4) + except: raise PybamError('\n\nCould not read from "' + str(self._file.name) + '"!\n') + else: raise PybamError('\n\nInput file was not a string or a file object. It was: "' + str(f) + '"\n') + + self.file_name = os.path.basename(os.path.realpath(self._file.name)) + self.file_directory = os.path.dirname(os.path.realpath(self._file.name)) + + if magic == 'BAM\1': + # The user has passed us already unzipped BAM data! Job done :) + data = 'BAM\1' + self._file.read(35536) + self.file_bytes_read += len(data) + self.file_decompressor = 'None' + while data: + yield data + data = self._file.read(35536) + self.file_bytes_read += len(data) + self._file.close() + DEVNULL.close() + raise StopIteration + + elif magic == "\x1f\x8b\x08\x04": # The user has passed us compressed gzip/bgzip data, which == typical for a BAM file + # use custom decompressor if provided: + if decompressor != False and decompressor != 'internal': + if type(f) == str: self._subprocess = subprocess.Popen( decompressor.replace('{}',f), shell=True, stdout=subprocess.PIPE, stderr=DEVNULL) + else: self._subprocess = subprocess.Popen('{ printf "'+magic+'"; cat; } | ' + decompressor, stdin=self._file, shell=True, stdout=subprocess.PIPE, stderr=DEVNULL) + self.file_decompressor = decompressor + data = self._subprocess.stdout.read(35536) + self.file_bytes_read += len(data) + while data: + yield data + data = self._subprocess.stdout.read(35536) + self.file_bytes_read += len(data) + self._file.close() + DEVNULL.close() + raise StopIteration + + # else look for pigz or gzip: + else: + try: + self._subprocess = subprocess.Popen(["pigz"],stdin=DEVNULL,stdout=DEVNULL,stderr=DEVNULL) + if self._subprocess.returncode == None: self._subprocess.kill() + use = 'pigz' + except OSError: + try: + self._subprocess = subprocess.Popen(["gzip"],stdin=DEVNULL,stdout=DEVNULL,stderr=DEVNULL) + if self._subprocess.returncode == None: self._subprocess.kill() + use = 'gzip' + except OSError: use = 'internal' + + if use != 'internal' and decompressor != 'internal': + if type(f) == str: self._subprocess = subprocess.Popen([ use , '--decompress','--stdout', f ], stdout=subprocess.PIPE, stderr=DEVNULL) + else: self._subprocess = subprocess.Popen('{ printf "'+magic+'"; cat; } | ' + use + ' --decompress --stdout', stdin=f, shell=True, stdout=subprocess.PIPE, stderr=DEVNULL) + time.sleep(1) + if self._subprocess.poll() == None: + data = self._subprocess.stdout.read(35536) + self.file_decompressor = use + self.file_bytes_read += len(data) + while data: + yield data + data = self._subprocess.stdout.read(35536) + self.file_bytes_read += len(data) + self._file.close() + DEVNULL.close() + raise StopIteration + + # Python's gzip module can't read from a stream that doesn't support seek(), and the zlib module cannot read the bgzip format without a lot of help: + self.file_decompressor = 'internal' + raw_data = magic + self._file.read(65536) + self.file_bytes_read = len(raw_data) + internal_cache = [] + blocks_left_to_grab = 50 + bs = 0 + checkpoint = 0 + decompress = zlib.decompress + while raw_data: + if len(raw_data) - bs < 35536: + raw_data = raw_data[bs:] + self._file.read(65536) + self.file_bytes_read += len(raw_data) - bs + bs = 0 + magic = raw_data[bs:bs+4] + if not magic: break # a child's heart + if magic != "\x1f\x8b\x08\x04": raise PybamError('\n\nThe input file != in a format I understand. First four bytes: ' + repr(magic) + '\n') + try: + more_bs = bs + unpack(" sam.\n\nThese two headers should always be the same, but apparently they are not:\nThe ASCII header looks like: ' + self.file_header + '\nWhile the binary header has the following chromosomes: ' + self.file_chromosomes + '\n') + + ## Variable parsing: + def new_entry(header_cache): + cache = header_cache # we keep a small cache of X bytes of decompressed BAM data, to smoothen out d==k access. + p = 0 # where the next alignment/entry starts in the cache + while True: + try: + while len(cache) < p + 4: cache = cache[p:] + next(self._generator); p = 0 # Grab enough bytes to parse blocksize + self.sam_block_size = unpack('> 4 , cigar_codes[cig & 0b1111] ) for cig in array('I', self.bam[self._end_of_qname : self._end_of_cigar ])] + @property + def sam_cigar_string(self): return ''.join( [ str(cig >> 4) + cigar_codes[cig & 0b1111] for cig in array('I', self.bam[self._end_of_qname : self._end_of_cigar ])]) + @property + def sam_seq(self): return ''.join( [ dna_codes[dna >> 4] + dna_codes[dna & 0b1111] for dna in array('B', self.bam[self._end_of_cigar : self._end_of_seq ])])[:self.sam_l_seq] # As DNA == 4 bits packed 2-per-byte, there might be a trailing '0000', so we can either + @property + def sam_qual(self): return ''.join( [ chr(ord(quality) + 33) for quality in self.bam[self._end_of_seq : self._end_of_qual ]]) + @property + def sam_tags_list(self): + result = [] + offset = self._end_of_qual + while offset != len(self.bam): + tag_name = self.bam[offset:offset+2] + tag_type = self.bam[offset+2] + if tag_type == 'Z': + offset_end = self.bam.index('\x00',offset+3)+1 + tag_data = self.bam[offset+3:offset_end-1] + elif tag_type in CtoPy: + offset_end = offset+3+py4py[tag_type] + tag_data = unpack(CtoPy[tag_type],self.bam[offset+3:offset_end])[0] + elif tag_type == 'B': + offset_end = offset+8+(unpack('"' +} + +wat = ''' +Main class: pybam.read +Github: http://github.com/JohnLonginotto/pybam + +[ Dynamic Parser Example ] + for alignment in pybam.read('/my/data.bam'): + print alignment.sam_seq + +[ Static Parser Example ] + for seq,mapq in pybam.read('/my/data.bam',['sam_seq','sam_mapq']): + print seq + print mapq + +[ Mixed Parser Example ] + my_bam = pybam.read('/my/data.bam',['sam_seq','sam_mapq']) + print my_bam._static_parser_code + for seq,mapq in my_bam: + if seq.startswith('ACGT') and mapq > 10: + print my_bam.sam + +[ Custom Decompressor (from file path) Example ] + my_bam = pybam.read('/my/data.bam.lzma',decompressor='lzma --decompress --stdout /my/data.bam.lzma') + +[ Custom Decompressor (from file object) Example ] + my_bam = pybam.read(sys.stdin,decompressor='lzma --decompress --stdout') # data given to lzma via stdin + +[ Force Internal bgzip Decompressor ] + my_bam = pybam.read('/my/data.bam',decompressor='internal') + +[ Parse Words (hah) ]''' +wat += '\n'+''.join([('\n===============================================================================================\n\n ' if code == 'file_alignments_read' or code == 'sam' else ' ')+(code+' ').ljust(25,'-')+description+'\n' for code,description in sorted(parse_codes.items())]) + '\n' + +class read: + ''' + [ Dynamic Parser Example ] + for alignment in pybam.read('/my/data.bam'): + print alignment.sam_seq + + [ Static Parser Example ] + for seq,mapq in pybam.read('/my/data.bam',['sam_seq','sam_mapq']): + print seq + print mapq + + [ Mixed Parser Example ] + my_bam = pybam.read('/my/data.bam',['sam_seq','sam_mapq']) + print my_bam._static_parser_code + for seq,mapq in my_bam: + if seq.startswith('ACGT') and mapq > 10: + print my_bam.sam + + [ Custom Decompressor (from file path) Example ] + my_bam = pybam.read('/my/data.bam.lzma',decompressor='lzma --decompress --stdout /my/data.bam.lzma') + + [ Custom Decompressor (from file object) Example ] + my_bam = pybam.read(sys.stdin,decompressor='lzma --decompress --stdout') # data given to lzma via stdin + + [ Force Internal bgzip Decompressor ] + my_bam = pybam.read('/my/data.bam',decompressor='internal') + + "print pybam.wat" in the python terminal to see the possible parsable values, + or v==it http://github.com/JohnLonginotto/pybam for the latest info. + ''' + + def __init__(self,f,fields=False,decompressor=False): + self.file_bytes_read = 0 + self.file_chromosomes = [] + self.file_alignments_read = 0 + self.file_chromosome_lengths = {} + + if fields != False: + if type(fields) != list or len(fields) == 0: + raise PybamError('\n\nFields for the static parser must be provided as a non-empty list. You gave a ' + str(type(fields)) + '\n') + else: + for field in fields: + if field.startswith('sam') or field.startswith('bam'): + if field not in parse_codes.keys(): + raise PybamError('\n\nStatic parser field "' + str(field) + '" from fields ' + str(fields) + ' is not known to this version of pybam!\nPrint "pybam.wat" to see available field names with explinations.\n') + else: + raise PybamError('\n\nStatic parser field "' + str(field) + '" from fields ' + str(fields) + ' does not start with "sam" or "bam" and thus is an avaliable field for the static parsing.\nPrint "pybam.wat" in interactive python to see available field names with explinations.\n') + + if decompressor: + if type(decompressor) == str: + if decompressor != 'internal' and '{}' not in decompressor: raise PybamError('\n\nWhen a custom decompressor is used and the input file is a string, the decompressor string must contain at least one occurence of "{}" to be substituted with a filepath by pybam.\n') + else: raise PybamError('\n\nUser-supplied decompressor must be a string that when run on the command line decompresses a named file (or stdin), to stdout:\ne.g. "lzma --decompress --stdout {}" if pybam is provided a path as input file, where {} is substituted for that path.\nor just "lzma --decompress --stdout" if pybam is provided a file object instead of a file path, as data from that file object will be piped via stdin to the decompression program.\n') + + ## First we make a generator that will return chunks of uncompressed data, regardless of how we choose to decompress: + def generator(): + DEVNULL = open(os.devnull, 'wb') + + # First we need to figure out what sort of file we have - whether it's gzip compressed, uncompressed, or something else entirely! + if type(f) == str: + try: self._file = open(f,'rb') + except: raise PybamError('\n\nCould not open "' + str(self._file.name) + '" for reading!\n') + try: magic = os.read(self._file.fileno(),4) + except: raise PybamError('\n\nCould not read from "' + str(self._file.name) + '"!\n') + elif type(f) == file: + self._file = f + try: magic = os.read(self._file.fileno(),4) + except: raise PybamError('\n\nCould not read from "' + str(self._file.name) + '"!\n') + else: raise PybamError('\n\nInput file was not a string or a file object. It was: "' + str(f) + '"\n') + + self.file_name = os.path.basename(os.path.realpath(self._file.name)) + self.file_directory = os.path.dirname(os.path.realpath(self._file.name)) + + if magic == 'BAM\1': + # The user has passed us already unzipped BAM data! Job done :) + data = 'BAM\1' + self._file.read(35536) + self.file_bytes_read += len(data) + self.file_decompressor = 'None' + while data: + yield data + data = self._file.read(35536) + self.file_bytes_read += len(data) + self._file.close() + DEVNULL.close() + raise StopIteration + + elif magic == "\x1f\x8b\x08\x04": # The user has passed us compressed gzip/bgzip data, which is typical for a BAM file + # use custom decompressor if provided: + if decompressor != False and decompressor != 'internal': + if type(f) == str: self._subprocess = subprocess.Popen( decompressor.replace('{}',f), shell=True, stdout=subprocess.PIPE, stderr=DEVNULL) + else: self._subprocess = subprocess.Popen('{ printf "'+magic+'"; cat; } | ' + decompressor, stdin=self._file, shell=True, stdout=subprocess.PIPE, stderr=DEVNULL) + self.file_decompressor = decompressor + data = self._subprocess.stdout.read(35536) + self.file_bytes_read += len(data) + while data: + yield data + data = self._subprocess.stdout.read(35536) + self.file_bytes_read += len(data) + self._file.close() + DEVNULL.close() + raise StopIteration + + # else look for pigz or gzip: + else: + try: + self._subprocess = subprocess.Popen(["pigz"],stdin=DEVNULL,stdout=DEVNULL,stderr=DEVNULL) + if self._subprocess.returncode == None: self._subprocess.kill() + use = 'pigz' + except OSError: + try: + self._subprocess = subprocess.Popen(["gzip"],stdin=DEVNULL,stdout=DEVNULL,stderr=DEVNULL) + if self._subprocess.returncode == None: self._subprocess.kill() + use = 'gzip' + except OSError: use = 'internal' + + if use != 'internal' and decompressor != 'internal': + if type(f) == str: self._subprocess = subprocess.Popen([ use , '--decompress','--stdout', f ], stdout=subprocess.PIPE, stderr=DEVNULL) + else: self._subprocess = subprocess.Popen('{ printf "'+magic+'"; cat; } | ' + use + ' --decompress --stdout', stdin=f, shell=True, stdout=subprocess.PIPE, stderr=DEVNULL) + time.sleep(1) + if self._subprocess.poll() == None: + data = self._subprocess.stdout.read(35536) + self.file_decompressor = use + self.file_bytes_read += len(data) + while data: + yield data + data = self._subprocess.stdout.read(35536) + self.file_bytes_read += len(data) + self._file.close() + DEVNULL.close() + raise StopIteration + + # Python's gzip module can't read from a stream that doesn't support seek(), and the zlib module cannot read the bgzip format without a lot of help: + self.file_decompressor = 'internal' + raw_data = magic + self._file.read(65536) + self.file_bytes_read = len(raw_data) + internal_cache = [] + blocks_left_to_grab = 50 + bs = 0 + checkpoint = 0 + decompress = zlib.decompress + while raw_data: + if len(raw_data) - bs < 35536: + raw_data = raw_data[bs:] + self._file.read(65536) + self.file_bytes_read += len(raw_data) - bs + bs = 0 + magic = raw_data[bs:bs+4] + if not magic: break # a child's heart + if magic != "\x1f\x8b\x08\x04": raise PybamError('\n\nThe input file is not in a format I understand. First four bytes: ' + repr(magic) + '\n') + try: + more_bs = bs + unpack(" sam.\n\nThese two headers should always be the same, but apparently they are not:\nThe ASCII header looks like: ' + self.file_header + '\nWhile the binary header has the following chromosomes: ' + self.file_chromosomes + '\n') + + ## Variable parsing: + def new_entry(header_cache): + cache = header_cache # we keep a small cache of X bytes of decompressed BAM data, to smoothen out disk access. + p = 0 # where the next alignment/entry starts in the cache + while True: + try: + while len(cache) < p + 4: cache = cache[p:] + next(self._generator); p = 0 # Grab enough bytes to parse blocksize + self.sam_block_size = unpack('> 4 , cigar_codes[cig & 0b1111] ) for cig in array('I', self.bam[self._end_of_qname : self._end_of_cigar ])] + @property + def sam_cigar_string(self): return ''.join( [ str(cig >> 4) + cigar_codes[cig & 0b1111] for cig in array('I', self.bam[self._end_of_qname : self._end_of_cigar ])]) + @property + def sam_seq(self): return ''.join( [ dna_codes[dna >> 4] + dna_codes[dna & 0b1111] for dna in array('B', self.bam[self._end_of_cigar : self._end_of_seq ])])[:self.sam_l_seq] # As DNA is 4 bits packed 2-per-byte, there might be a trailing '0000', so we can either + @property + def sam_qual(self): return ''.join( [ chr(ord(quality) + 33) for quality in self.bam[self._end_of_seq : self._end_of_qual ]]) + @property + def sam_tags_list(self): + result = [] + offset = self._end_of_qual + while offset != len(self.bam): + tag_name = self.bam[offset:offset+2] + tag_type = self.bam[offset+2] + if tag_type == 'Z': + offset_end = self.bam.index('\x00',offset+3)+1 + tag_data = self.bam[offset+3:offset_end-1] + elif tag_type in CtoPy: + offset_end = offset+3+py4py[tag_type] + tag_data = unpack(CtoPy[tag_type],self.bam[offset+3:offset_end])[0] + elif tag_type == 'B': + offset_end = offset+8+(unpack('= 2: + #print(allele_counts_list) + if len(set(snp_genotypes)) == 1 or allele_counts_list[0] == allele_counts_list[1]: + # If only heterozygous sites 0/1; skip the site (equivalent to n bin or n/2 bin for folded) + # skip if all individuals have the same genotype + line = inputgz.readline() + continue + if len(ALT) >= 2: #pass count_pluriall +=1 # TODO - work in progress @@ -95,10 +114,19 @@ def sfs_from_vcf(n, vcf_file, folded = True, diploid = True, phased = False, ver # SFS_values[min(allele_counts_list)-1] += 1/len(ALT) # allele_counts_list.remove(min(allele_counts_list)) else: - SFS_values[min(allele_counts_list)-1] += 1 + if folded: + SFS_values[min(allele_counts_list)-SFS_dim[0]] += 1 + else : + # if unfolded, count the Ones (ALT allele) + #print(snp_genotypes, snp_genotypes.count(1)) + SFS_values[snp_genotypes.count(1)-SFS_dim[0]] += 1 + # all the parsing is done, change line line = inputgz.readline() if verbose: print("SFS=", SFS_values) + if strip: + del SFS_values[0] + del SFS_values[n] print("Pluriallelic sites =", count_pluriall) return SFS_values, count_pluriall @@ -163,20 +191,39 @@ def sfs_from_parsed_vcf(n, vcf_dict, folded = True, diploid = True, phased = Fal return SFS_values, count_pluriall -def barplot_sfs(sfs, folded=True, title = "Barplot"): +def barplot_sfs(sfs, xlab, ylab, folded=True, title = "Barplot", transformed = False): sfs_val = [] n = len(sfs.values()) for k in range(1, n): ksi = list(sfs.values())[k-1] # k+1 because k starts from 0 - if folded: - sfs_val.append(ksi * k * (n - k)) + # if folded: + # # ?check if 2*n or not? + # sfs_val.append(ksi * k * (2*n - k)) + # else: + # if transformed: + # sfs_val.append(ksi * k) + # else: + # sfs_val.append(ksi) + if transformed: + if folded: + sfs_val.append(ksi * k * (2*n - k)) + else: + sfs_val.append(ksi * k) else: - sfs_val.append(ksi * k) + sfs_val.append(ksi) + #terminal case, same for folded or unfolded - sfs_val.append(list(sfs.values())[n-1] * n) + if transformed: + sfs_val.append(list(sfs.values())[n-1] * n) + else: + sfs_val.append(list(sfs.values())[n-1]) #build the plot title = title+" [folded="+str(folded)+"]" + if ylab: + plt.ylabel(ylab) + if xlab: + plt.xlabel(xlab) plt.title(title) plt.bar([i+1 for i in sfs.keys()], sfs_val) plt.show() @@ -189,5 +236,5 @@ if __name__ == "__main__": # PARAM : Nb of indiv n = int(sys.argv[2]) - sfs = sfs_from_vcf(n, sys.argv[1], folded = True, diploid = True, phased = False) + sfs = sfs_from_vcf(n, sys.argv[1], folded = True, diploid = True, phased = False, strip = True) print(sfs) diff --git a/stats_sfs.py b/stats_sfs.py index 4c376fb..b1c5128 100644 --- a/stats_sfs.py +++ b/stats_sfs.py @@ -1,4 +1,3 @@ -from frst import vcf_to_sfs import math diff --git a/vcf_utils.py b/vcf_utils.py index 049a6fa..7a532b6 100755 --- a/vcf_utils.py +++ b/vcf_utils.py @@ -214,6 +214,25 @@ def free(obj): del obj gc.collect() +def random_indiv_from_vcfs(vcf_files_list): + """ + Either for two VCFs of two indiv or two samplenames from the same VCF. + """ + if len(vcf_files_list)>1: + # two invid from two VCFs + for vcf in vcf_files_list: + with open(vcf) as vcf_stream: + for line in vcf_stream: + if line.startswith('#'): + continue + genotype = line.split("\t")[-1].strip() + if genotype.split(":")[1] != '0': + genotype_list = genotype.split(':')[0].split("/") + print(genotype_list) + else: + # two samplename to randomly pick + pass + if __name__ == "__main__": # check args if len(sys.argv) !=2: