From 4623c0ea8979b23c82ac8493ffabb6d35594d8ba Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 27 Aug 2026 09:17:00 +0200 Subject: [PATCH] test(determinism): lock output against map iteration order, and make the SJ sorts total Junction counts live in a DashMap, whose iteration order varies with hashing and with concurrent insertion. Every path that emits an order already sorted, so no order was escaping, but two of those sorts keyed on (chr, start, end) while the key also carries strand and motif: a tie would have fallen back to the map's order. They now sort on the whole key. tests/determinism.rs is the lock: the same reads at one thread and at eight, and two runs at eight threads, must produce byte-identical SJ.out.tab and Aligned.out.sam, in single-pass and in two-pass mode, including pass 1's own SJ.out.tab. The test was checked against a positive control before being trusted: with the SJ writer ordering rows by a per-process random hash, all three tests fail. A deterministic permutation does not fail them, which is correct, and is why the control had to be per-process random rather than a fixed swap. Answers #210. --- .claude/worktrees/agent-aa8bc29ee5b6ad4d5 | 1 + CHANGELOG.md | 9 +- src/junction/sj_output.rs | 19 +- src/ruSTAR.code-workspace | 7 + test/__pycache__/bench_report.cpython-312.pyc | Bin 0 -> 4072 bytes test/__pycache__/nfcore_diff.cpython-312.pyc | Bin 0 -> 9827 bytes test/__pycache__/speed_bench.cpython-312.pyc | Bin 0 -> 6031 bytes tests/determinism.rs | 276 ++++++++++++++++++ 8 files changed, 307 insertions(+), 5 deletions(-) create mode 160000 .claude/worktrees/agent-aa8bc29ee5b6ad4d5 create mode 100644 src/ruSTAR.code-workspace create mode 100644 test/__pycache__/bench_report.cpython-312.pyc create mode 100644 test/__pycache__/nfcore_diff.cpython-312.pyc create mode 100644 test/__pycache__/speed_bench.cpython-312.pyc create mode 100644 tests/determinism.rs diff --git a/.claude/worktrees/agent-aa8bc29ee5b6ad4d5 b/.claude/worktrees/agent-aa8bc29ee5b6ad4d5 new file mode 160000 index 0000000..40d9286 --- /dev/null +++ b/.claude/worktrees/agent-aa8bc29ee5b6ad4d5 @@ -0,0 +1 @@ +Subproject commit 40d92861cf3a1c2125d68a589655c0277fa519cb diff --git a/CHANGELOG.md b/CHANGELOG.md index 7c20e6c..1a64eaa 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -20,7 +20,14 @@ Sections commonly used: Features, Bug fixes, Other changes. implementation built the SA in 172.953 s versus 267.592 s for its original 0.7 baseline: 35.4% faster, with peak RSS reduced from 10,512,408 to 9,169,892 KiB. The complete output hash was unchanged. - +- The splice-junction sorts that produce `SJ.out.tab`, the `SJ` solo-feature + rows and the `BySJout` survivor set now order on the whole key (chromosome, + start, end, strand, motif) rather than on coordinates alone. The counts come + from a `DashMap`, whose iteration order varies with hashing and with + concurrent insertion, so a tie left to that order would have been a file + that differs between runs or thread counts. `tests/determinism.rs` locks it: + the same reads at one and at eight threads, and two runs at eight threads, + produce byte-identical output in single-pass and two-pass mode. Answers #210. - `cluster_seeds` reuses its window-bin map across reads on a thread instead of rebuilding it per read. Merging two windows re-keys every bin in the merged span, so the per-read pre-sizing was only a floor and the map diff --git a/src/junction/sj_output.rs b/src/junction/sj_output.rs index 76b6527..05c3082 100644 --- a/src/junction/sj_output.rs +++ b/src/junction/sj_output.rs @@ -165,12 +165,17 @@ impl SpliceJunctionStats { }) .collect(); - // Sort by chromosome, start, end (for distance calculation) + // Sort by chromosome, start, end (for distance calculation). Strand and + // motif join the key because the source is a `DashMap`, whose iteration + // order is not stable across runs or thread counts: a tie left to that + // order would carry it forward (#210). junctions.sort_by(|a, b| { a.0.chr_idx .cmp(&b.0.chr_idx) .then(a.0.intron_start.cmp(&b.0.intron_start)) .then(a.0.intron_end.cmp(&b.0.intron_end)) + .then(a.0.strand.cmp(&b.0.strand)) + .then(a.0.motif.cmp(&b.0.motif)) }); let overhang_min = ¶ms.out_sj_filter_overhang_min; @@ -274,17 +279,18 @@ impl SpliceJunctionStats { /// keys so the SJ recorder can be mapped to matrix rows. pub(crate) fn sj_feature_order(&self, params: &Parameters) -> Vec<(u64, u64)> { let surviving = self.compute_surviving_junctions(params); - let mut keys: Vec<(usize, u64, u64)> = self + let mut keys: Vec<(usize, u64, u64, u8, u8)> = self .junctions .iter() .filter(|e| surviving.contains(e.key())) .map(|e| { let k = e.key(); - (k.chr_idx, k.intron_start, k.intron_end) + (k.chr_idx, k.intron_start, k.intron_end, k.strand, k.motif) }) .collect(); + // Total key again: these are the SJ matrix row positions (#210). keys.sort_unstable(); - keys.into_iter().map(|(_, s, e)| (s, e)).collect() + keys.into_iter().map(|(_, s, e, _, _)| (s, e)).collect() } /// Write the 9-column `SJ.out.tab` lines (sorted) to `writer`; returns the @@ -315,11 +321,16 @@ impl SpliceJunctionStats { }) .collect(); + // Total order, for the same reason as `compute_surviving_junctions`: + // these rows are the bytes of SJ.out.tab, so a tie broken by `DashMap` + // iteration order would be a file that differs between runs (#210). output_junctions.sort_by(|a, b| { a.0.chr_idx .cmp(&b.0.chr_idx) .then(a.0.intron_start.cmp(&b.0.intron_start)) .then(a.0.intron_end.cmp(&b.0.intron_end)) + .then(a.0.strand.cmp(&b.0.strand)) + .then(a.0.motif.cmp(&b.0.motif)) }); let mut written = 0u32; diff --git a/src/ruSTAR.code-workspace b/src/ruSTAR.code-workspace new file mode 100644 index 0000000..9e68e72 --- /dev/null +++ b/src/ruSTAR.code-workspace @@ -0,0 +1,7 @@ +{ + "folders": [ + { + "path": ".." + } + ] +} \ No newline at end of file diff --git a/test/__pycache__/bench_report.cpython-312.pyc b/test/__pycache__/bench_report.cpython-312.pyc new file mode 100644 index 0000000000000000000000000000000000000000..d9ad800bd3a71e912209a3e4e1fa58dae6f4578d GIT binary patch literal 4072 zcmbUkZERE5^}hG6KVmz!^VyU_`%3tTAxSB;e5Hh@38jS=XrWsd!mMZCm&A$vZ0>yq zVm)uBn^dRGgtOA&mZD{XA7QPu{6S*UexUu>pQ)08v@@@nw9!edKc=}#t4Y&-?A+(r zF^rCGSMt5*oqO&%=bn4M{QLTP7lM&|D>vqFBlJ1{s0>>f5EhF7oI?Z+A|enmD#XM= zQQ)u?l?Eki#e*1D9JR&lgLVOlNR`IC6-LV|v3)3%zXqM(NMo$9T|98c02?9@0p>u-2>ZkKC-^N|x27%<72sox6QDb8JJ$(Y~%hr7N|mYe%sA*v7yR+jw+{05rrlbnM=# zG2^#am=R+JgRKp47Ul;Mbd)LhiH=$h$sG^Wp{zn#&OP&;cRyYCS9f9=U>}v@=KyOowF;igSCGWLpFwC$mES>=DITS%@ z!9ZiR6c!mm70+X$AsmLdNsBolc&atKFt0k%M89>M23DS~ouk4Q1RiNe7 zsPal<1JgzWk3odrI|!~wOp8;Ux~H=L0FE}p#)PQ>!#3utoLF8fdJe%I)AqC@?KDK{ zF+{^Yz73+(VTk7uu^IOB=z^pPhCN-UH&jj=j*32dV-@Fo7{C-*#1M?SX8@i$jR8K4 z0`>zYG96S+M3t}#4r*5^jZ^joe+5%`BH#kdsR+}VSw~fc1aOzG8w-!<(P3>3rBxT2f3!~K0P?E!{30Wn}okR)gumx4JX zQy4LYn2C8DoAz=Xn>eXObW<8vCm7GNa&`Q(Sy!3lA$eABZkE8KE>ip?$O_D=83_r0?(yQk=BUhw$Nyg2z{j=X&k_I-Xp z6<0m_^Vi>d{ruo{=ZY-O4(G_N_Kw-t=EyIP}=FM6&wb!Vk)|69&d3u@>swIYu%vu_a;NB57FmtY5Ki2y$Q5JE#Dc!(7bhB!V1 zR~Du`1~SAHcq$RkV_mVKwo5Kj+91+a@z}mB25G@Y~-~hN6@TE?e0OoE^9c7ps zp7J;WJP#7e<0>U}Zyq)Uc?4jje*AG?N|H5p83iQM!PIbEBaE(rEz`!rN>pV$R^jpJ zI*yZuV>;~y=n43-y)b3aXAd2ibc5a0qPsDt-gG}wY+W^Dn~u+f3oRRFAD-&%;@s0CzN=FzM9 zj8t<(e;bx=Y4FecWk5r$ z*o^1i%*q-GkP9LE)n>e^?tph)sL^<88{}=LAsILYZVWTrR}+XkZPS~ox>f0k$Fj3X z>Vd|awts7xwYpTR0oT#65#M>t^EB~aa8+(3ps63>TAMX2@6^a_T*lt@B2CNq8tUd6 z?EUR7bGz4^_)j}7y-9_WrFGRhp!I4Dz7c3SKzXtGf=5v5`S$j9`8l43p|cxjF4uv$ zY-MHsS#Hg&oJ@LTrt;27m6Hi52Sg>%g`HgL6uYU#tnF@Qa>r6T=Sp=NHCzo~1Faxc zhh;!7b&UXqp`9C1Wh2#Kbai#z^-GC{3@OTxlf0*f^>*eg%$j;qHmn>5STyqfdwl$Q zcan#@eEa>Gzh6B=`GFMP*!|?%lnA_^b@f#QDIw?zW!?vZHoR*bQQ?|#O;(`YR4Bh; zQqy>!36&GNVQ6xpuS{)7(0DAN%g`HBwF^k@YDYPVd_tup64tv@t#{wMY!A?)u6jGD zr&3haOo>e}6DxG&gjN11&uUN)c?)gY8PrFWn9_(^PQ0e1gjm2~;!!B+CgwJztN8(8 z)D%Ldji^ywF|pMNTec_{-g5fJG3eN7FOYAAA3F&-HG>w`Z}@1-ueM*_elb0F{Koq2 zQ%|35J{vwaIz5^{K0SW3VfDPb{c0#P@NMb-?Dj1<>kHo3=bhgx`uu;vbp_YjdA#nH z!)@W8nRo88fXC;Zs~4JAYo)ZK0mX4wztr_bKc)m^tat=Uz2+!-!T1JUis<%+4k(71-X4@&qdelksI=p z(>vy9b^w}!+;cfSH=OUkB+c|)wC7LFu#dLP9{+IbY{y*toO*R&Zs^*YYpvI4p>f|Q zp8cf;)bjXOji}CB@V+wdeAQBS=e)D0j<8Dy{&?W4< z=qo(7X-+IOZu!LXR7nEPe{y@QdS^O=LH;7aU;xuq@FL}>w^LPzzho&tcci0Cr=8%U z(1slYs!mSz2%nyX30y~~P~6Kf_)EZVI)XtG4+n#$D;OM30wXmTq)bg=758D zanVSqEY1-IpN>RT$_ox{uTPB-!$UbV}v<$NXLume7mt7Pfoa zgC@19k$2ca6uUtaW@@zAO;*WJDL(eYseHuy_@FBJkjI{&X`ISzc7FIX%cf>)@{@C0 zYDs22o2^u(jqdw#?!D)I_!q0ygu&OhkO?oeVAyX^Mf>!n0{-GlDBQ$kF-A?3CYVD_npr((et@xtkBD-KhLVvr zvE~oSNeheqtng=J?H^K;R<@CKe3zWGv0GRt*8ZRfxE$o9Q3n35h`M*$uz)3Lo6HNI)sqK zu}mN|D<*lapLwO5>F(}k621@*%{=F0h05yt#yc55?0YjLgnZmS#y`!UVG@y~U_$$y z9zM^Aye}pMxfpkVnVAmxr&8GC0kMaVEft(BX@*Oen^3vyd`LDDaFA zK#C3RK0e6R%!qf;C&sO4lVB&r7!7{4zTi;DRlITZ={nM)VX39m8$hc@UY7h~ZKSeVVgmw@T@ z$IHT7b|&VSjt)L4h(4a_?1D7NiCsJw0XX3S)%FfwBovInHC5 zxKkzw*9R|i>WnEIpXy}Z1~9oqoEJOdvB=!pj4v|dn-drgc0wDViwpy}Oo9)^paos3 zKORl^_z+xp6)8+J(>{2)IFl4OgeeiU2f7J_)qV~^JH@FN9lHs58$wXPB|by6;K3EM z`q(}(6psm`P$$p(#A#*kp@WC~M*_a7L%m144;=*5<@%<2x)1jq=8pCp?d|LFPX&%n z_4W9qBh#Xo5c<2ic;8HCFeFYVrvMP%ACHMZRGk1Z+-He;5ZQE;$wugs8sSySKz0G1 z$D>?lzz0_>?cUhI3Eml4*ULo{;+$6uv7JFN&>57Z;Qw{&5*c?Pcn~;Ur^knSdVo~> zy2dYFJk_@R}D(X&ZP0K0=KioM95?o}#PvZCFZk`F zqE+i*Ng02DeWU}(sT#;_T7TG~m8yBzRV`~JStk?K`M1L)#u6+RG0l@|Zl0%DJ)D$- zlj&s;p|TzXDTq(?tOnWeN6(7iZk5S!leSkiudHM&t;Cu%7+%5VX|XaxHn8TJ9-~Z! zo2#p6C7G6u+9(kC>fbllZwE2i%JafC!!T=$;{5dxrWIeTix$eTmE@V};xl z9Y3Y$W~Mwgg^FHhL%c!_pE)-)u23P7iwcTf-AGVK#P<~g?1==tLa7W{(I@#xBs8TM zcn)Yv5EU~t@SMnpKq@F?Pzog!y)c~=LlK3HC%Bkm_Qw-*0btPasj$bSkVxVvWXu=k z6awZdBn!*%+W}~vL6GAy*ee3sh(ZJ)Ap`NGC?F(70;}|W?LtiKCHV)pfb;_X01WMh zpYY$1q_84MzjyZCvl%hRrp~UC``7BX-oJA9%3oevFs8}0zi4-5*oE(;$)eeQ=hcN* zZy(Q2e0=$X%MY$RY3*Ka?S9nrX=`86;$EYz8Dn;5wxM8S^5kc9#~N+BV_UF2q1%?} zwp_2;wr<86+F!Ua%JGzR{>JJ`<3)!%+mOSvLT-1V(X;IEq^Y9aeMer9zyD6rS%1Ij zZd0}|H=OUux92Yx+72%_A6|5?I!`RM|4Mv3R2YAA_4HeXt8W(_4Nn|x%MNv;RY%8H z1U71<;-N(cV_mY8UX*dovqZ5J4l|CThKlp$<>1w!IS(P~ zu|j*_qf?8-Khyt6FS-87`smhD;PE(MtMlS#^!Pe~(Jfz5uta*>*8+N~@9#ND62Bm+ zlU9Sj7FCewP(|h6s;DO8vJM^q|465@dJ-6U1^HlP`D1ktiR$B2k+3|QEX~$N%u_NJ z#6ge{0xjzi0t9suh_rcq3p+dFIp>S+~!dVN&hMv0G9@fRW*?P7C?8_|( z&Y)F=&`<**pu2kOO)o+lmg zfi6I>nQWO;nOCx8re)dGnmSb^4ey$B z*>)3qXHM2dVfB#s4Y=t`Z5*(Wa#8KCkXp5$tuP?i?KLANRy26I z(gUPR{IRu;#3&7wJTj43~MN z&@1>ZkG>?PC?OF6=|cf}gN{doUqV0zDxh&pR7`L9B1vwD=i|JG;?JR(Bp;svkrfHW zKzWS#rnrbg-9)1VbRRY7ZB%#l0txZgpeG7HVF40A@#fal;2Nk<`pj&0zEI!sZ;no+ zso6kfdkvJbVq^1t>8_M(&)ryQ^ekA?R61F7)ZaJXHQ%#>re${AIk9je+rDbv2BX(? zgweUywkyZ|cwc%jQ-8ntZu7m?Pc7TmotU$0y#cd0Qs4gl?~kN=uPQA0iTfmBtRtA9WKc_}jYcKddt4^Lb%|}o0Nz(st4zNs z?whm~1%S6^KRpOZn634K*u%p}EqLhz(Pz;#K$b!R74X}Da2RZ*AClPZ`y^&m^~73D zv4qx{-U_xwLEMdK#SLuMCK!e_dL;(&if=XAyedP7y;kGb-67}f=Ar5?GV9Z2( zLH-PC(#57XVf7YpDpWj%JR*g8_zC|CiAD}aYo;MXXYB>s-c_R~f8nXSK0OF#%$-{c zx3a>jy}jtzQmk)6Mo!MXQos92{ek8B1Nn)S`h)4=HCIFCX4ZFi?%oM?Qil7UoE3B2 zU&(p#k^k={G`F!eJC*g{`%Z=^x*G1ExO*bE=QG#7wWgNrc+Q!t%f9)bC3CX4z5V0a z4`%ba{J?`-kDSZf4`v3lb$3TTcQ-vZVq4mu+p+ovq)F~u->MSUSIyOOiK6Zfbe%j* z`~u&P(!8hR^S*Nvvr>SE9 z02OZHv5|EsDDlwo38Du0=7|k65F>o)gBeb?u4V=)A~i_;z|&LaPc>~7`8BVfr^Loh z##Y;^V`<h1#)6kdvLQgdfH--_~&XbmEG z5fgBb>8+3z8)Lfti6w{-_#paGh*H`Zx>rDXR^x<3FNM}r>vR}^G?|*n`JBF zZB@u+8{3HBl_C72=8kgiWLp4kgSH01?G(3FcCexwLEbvxx^CP1MBiuJ0-mSr2&H3XSXmOJf+6rM$b(U5V_h_Zf z!gqij6d0sk9};Z$!yXO$wDkZY3*TME9QI&s58yP=^?ZFzf3IARJh#c1|eZ6ZU0|lP2ij~uw1|gww0tWg_aZ&#g*OzP4KA$1#n6hTqQyT?hz41xkNt6!D@P+9|djD^e_a3r0p|2xGhVuH;8pX1QLZJ zgOCdh1f$Rta+-@IRCgwi%+|ATT5yl5@DEG%qiQsei3cEzm`JMeKUVON9vzPuvD9*@ z96V6Nft7$}Oq~UTklrWlsV<_Lo7z|iPpMIc1VpDIaeoNLBRyVfe|fwXc<6_SR-93z z4=@&K?9%pc9;*f<__tw{qK_gMtY0FI^>y<-XdZE*r(2;|E&|Vl#GXW*INI5}OCtIi zSt6JYMj~aVi-F>CNYozY1XEJ_L8TAysw%XF<+ud+S~mnX_l}Nu)a4Ea6ns{)Xk$Qa zVkAc;*o&TdkHAPeOZDosNIhWV7Bm?lOp4vmkDof{knG=wm4 z4BT(XSpe?$IJnA5gpi6!MM$s6&4R7OCZh>KF_!|KUKH^3G$=YBV515MVycQ65t6s; z2C$a!Qmf+N1>pg8A=U);6iN*-@OuEB@?iuVQj--L`4v`gV8`1x=0 zGpIrd5at3f`ErJWR;qF;G(fj%iTy!R@;s1H@Wb;Tz^FU$6aE7b3AqrRPmIj6ky$ae zrA`(}{d;5Yj%Cgk_MS0bSmGu zh%fG0YF}(Em`<;fBTpMz7EIv#xKwo3-LJn}pGjtY_qOMD{B+lkcIB!3jURWxI=VH7 zD{~UmZm?0ew5I8dYr%>F4E5kRF=u3Pz`Z(CMKj_c;@4mKTZUdY4&e+0OHkl3P#`AS~GJook4)plGMZ=@HC1Pplar@G% zkNv9`uNJPoy?oJEnBoc-0)>{~XO`)8E7sWhI|t@&gviK+Ro9+}uKdJe&#LEm`jyns z=dEq&q0G?tUVljCxlb&qq2hM9k@TM}gJX#_7OZ>nM8R?(Kal6(Ca%Cu>{=ogo0kTb zxW_|F*9sTM3#N%x@>0=Yxh-Wom+hWU4W1(HxMN?iXAk6(xnQB=_#(IX=F-WfxB!KR^NEL;EON6ktkfcQMe%#wuq~A5_zHL{(NqYwxlQSTwb{RzP%K0`AOeT`+wBG zc;xYg!sM0Z&Z`A3u-YD6rKf*wu){^BudmWupE{e;18Wv%=49qZ#-G{y{V{M&*KPg6 zz2%8}=cn$S=^F)le|`rf@W=%dvr6tLnl1c!L@~4%NF@~WdcmoOh*V`6D5i$|@d`N_ zk8xlW9)%vs7y(o73RrvB_fd4N6_@x6eJmPhlM(I&|2NQuYy|qV_HS8z+~I;`oWu9QAEQbt{h z0{2d9Td`?dv6(^t+>Ou8dYdj~dQs=XJ6{|i@d^Betrg$$;#;^Ccf4TiknB8+Ti3xW iNLtd;r$kecw4_eGH}dYt!o>S;d`7wfoye`~Vg3i%g}2rK literal 0 HcmV?d00001 diff --git a/test/__pycache__/speed_bench.cpython-312.pyc b/test/__pycache__/speed_bench.cpython-312.pyc new file mode 100644 index 0000000000000000000000000000000000000000..ca6d35d2a58e5f7669f5352003844717baaf4004 GIT binary patch literal 6031 zcmcIoU2GfIm7d`Y$>ER`B~gxS$C9ohTQ(b!vSK;@ch;8VpJXE&mYmqCQ+mWZq)3r7 z^vuwbW{0s{1Td8ZFths*6;4nEYal8TVD+#MDt3XYzydY8FJz=b>_lA{XkPl}z$ps! zspk$kqU@?|UUqEFJ@?#m&pr3dJ?EZ}|Lgbj6g-~yi>bawiuyh2*bi4F=%?K@MSV%x$&Y7zZtz_GI`mo~A^1>IHk5 zdYjr+r)>(;4odXHA1^k10d|ll)lqG?zQ3q-QfMo;O55|ka8&!2?ap*lrg^LxXb*P1 z$BqRtkKNF=h&wPhd4Qsnov^Aq)$`Z-fZt-TXX~b^8D@qFd4&&Pk++3(8fnRlhy)3v ztSC&PjF^#?85EZlq~ughP{KkwnUF*U2?-%7sTvv?J#(c8i7-KvNfY7utE#9XZ9+t< z0J@}jI3_D%uObPmcnigo(^^h}B~)}%ij$EjPR8Sr?95Gc;P63o@Zdp|6_SdGks=CM z?cqlvC?+c6?WCF%L=>A)E~0EYS6$)!xywjX1W8SZk{IeC)+eO|N=jIqCIs1tWpPH7 zG+q_6f+A>QdIpV)F(Idl#LO8afIY+rSVhv%t(>T8Nm){l*YzY-gk?|-lQ9~f;niFm zHWC&TMVLV`kkyo2j7;jOF9&i}FNiRQ7y??g*aS+;a#jU5A$YPg+Jr0}LK-aDt7btm zIxb4F32+ArhpR3}kzR0oLezQ{F)e}Z)Y$R)sm5i}rY z;6sQ+Rz)&cYoC>j5SvI!Vub$#Ildv<^a>J86s*rN41M;ppiMyQjEmsPSS>QFvrgik zp(y(4EO0-jv|7|b+-;SzL@7$EwJ^$VlADAb9Vyo)Ym>-3G>_d$Sy{;LRa-m4@{H!U zM+_R%b`&3`@~jrHM^cczo_`F=V6_IDn`*YHdi$m03)m@n@=nZLTBykEQQ1&+i8P1>| zP@lSN4)~;PQ_gKTY+vO$Z8Mvv+HHw?i#?XDPJ_d|?xiS+1Z1B6 zm7aooz(Xq3FPci~l@$ zWoX#+D7k1xm{wHPWXW}4dNMLZp)4h1CM#z}$%Onjo>k zjmcOv8FdB}kd$GvDM$h)14>OUn-(;>vYg`1LXgbLqqnUKvWQl&0 zCGIkvLR=FSV8NASy2(458COa0ksMlm{v6yUhnh}eKUVgGDEY$bxA2;!zVYu^_ji>2 z9S^y6bhM0)uKJIy`v=PYf!T8#zQ+5*bHl5?J+ng_9)BTLJoI7veBVZ3*Zp_r-kpyv zcK+w<>-|IJ{-L$NxdQtqhT?aA;|<(DJ$HK5`{KM<_IAw#!qf6Ui+%;Z~+%b-=AFWUpY}e zacTB^#TzKL7EYDeoqq~~ZH4wft3KHEFB|p`_A~#}?;JekjzNF(A-TxPnNr$uV(|W$^E?D6t%fCx}HG|dfNU|pY0VL$eN8L%@OWV{{rNNs*&b-rz7o->@d2{+^57D=}~id(jD z9e4*>d2{Rq+=_R?4mNYAy%I&?U63Jn>%d6(^B!%pE;c+hD_*QwQM-|EK%3|DUc;H7 z^FG6+xD6jzwI|`w!Lz{e{SEutp0Tfet9{RmbZi~L2zUNH%M;$}s$1;#)$D$08@v7a z27@=ehTmwwU7xT9W3VS-=N~?1DQ)wH!h3DM!8fUM_igV#_|Dy*GTYABzb*Gbor^Z?m{vE`LI7yO7| zlg~ZMv1ei!JvR(pTSKbrugA1bo0r;aQ}q`7dg#l3d(XBWz)^=ouMX>+U9{^R|DbAz zBG!$?xR6V0x+@$e_nyfTFhzP(vNEYR^=g@{UF%!eK=*{hq&f#!WYYaaFr~`UInz}I z1fwP!P2y=N8Qi4221sCXaQ&)URGu{HD47j+^_%2wZUDR_XTng}4$z2R_`GniEL zX1jAv0JF|Pd7To1ZSo#iIV&sekwT8z5mViyibf2}xsAc4^9t?+}av>Uq0SaDL zL;y`D3pPe|K5EU3klUZPAlenOFWnCcfUs7*DCvG{x6wXp;ahh1m0A~81N|ds-qdEY zP;dYNdM#x{Z-_Ph)*{hby=8Z66@=-xMC+r~Zf#U=?ds}6LlO>aa@fMR7PeG*z3*HV z!ybPg+mZ?d5zy0vwjfz3^K>3Cj%!JzO$h6SFXso^{z2+~1@-eiQs;b;Cr zj*vrn6^@}vt3VwEN;yutOj3&L4wy>HA)NtjIuq%O!|N4zO@p$5%LK_=MWzaXUQhvk z10anlvRVbvQvmjAB`2y?bts|8xvT{slK>P|{emn3YKGaQijC*eXlg>N0XMQ2Kzmg@ z2XGFsJ%M_p@&Yt$Em4Ko6M(NHe;~OZ4!?_h8A-E>*-+5r0O7)+RpDr!Q4%=<*DsT< zV)_J(qcCdmOx_YAeN`dI4Jt&46@`={CQE`yA-{Z>&delEDo}!Iis`Po-l_skXO@6v z0FY@}&`h@?LOlxgFiRrI1?o`!ru4Zc>Ho03R2X)$o_}tXP%|QrH*6E7gq4ECVx3u8hN{P;@WEa^)k1+ zz!XM`2PzHC^Gq=}KUlm`Zs;s96)289bDsIVi!BAu8W-7UYAKG+tMlSQ-~7e$j(w|5 z-38aC~1K=&e33AD{?KQ(@0R04Y! z;)}OGlb?87_cx=(4$=G>a|Cc zrRz6J4e$Jh`_U6Owete~9S7FFciu7u&eUa#7+qtpJ_a;2^W&Kf&Q}<{e|_%y2Z43& zr84)@!~V~Ye|CH^u`H~0{nHwErsDNIg(FRwTokZrQi?`HE@clebWa*?xF15GzoQVC zU(E`xDh}0yS}vOwe`p1J53-tfAbpM#O^`XYP^;;VMzI`=Mom5%jpxX3lhLRG4YivTEvl-&yX-9jP7JGo=EiF9(jmWD`=BsEA`nleNLz|uI}e4v@63^{#Ugv?;w z%cghBFP~QKvT~CZxMZ3Hmn9fMei`{0rQC$38s6MVNHV#!cv_KxM_`EhFYtQ8&@}xm z>!aB}?xE(n;c{zu zwY7J>_0@9gtE<5yYt+%-QEg9MetOqabd2^tnWib$SJ1y=c2!v4?7+PXcQ4G1e(=_B N*p?>@#qPF@`wL~g`Jw;- literal 0 HcmV?d00001 diff --git a/tests/determinism.rs b/tests/determinism.rs new file mode 100644 index 0000000..a778748 --- /dev/null +++ b/tests/determinism.rs @@ -0,0 +1,276 @@ +//! Output must not depend on how many threads produced it. +//! +//! Junction counts live in a `DashMap`, whose iteration order varies with +//! hashing and with concurrent insertion. Every path that emits an order sorts +//! first, and these tests are what keeps that true: they align the same reads +//! at one thread and at eight and compare the output files byte for byte +//! (issue #210). + +use assert_cmd::cargo::cargo_bin_cmd; +use std::fs; +use std::io::Write; +use std::path::{Path, PathBuf}; +use tempfile::TempDir; + +/// Two exons with a GT-AG intron between them, so reads spanning the junction +/// produce SJ.out.tab rows rather than an empty file. +const EXON1: (usize, usize) = (2_000, 2_300); +const INTRON: (usize, usize) = (2_300, 2_800); +const EXON2: (usize, usize) = (2_800, 3_100); +/// A second junction close to the first, to give the neighbour-distance filter +/// something to compute and the sort something to order. +const EXON3: (usize, usize) = (6_000, 6_200); +const INTRON2: (usize, usize) = (6_200, 6_600); +const EXON4: (usize, usize) = (6_600, 6_900); + +fn lcg_seq(seed: u32, length: usize) -> Vec { + let bases: [u8; 4] = *b"ACGT"; + let mut state = seed; + let mut seq = Vec::with_capacity(length); + for _ in 0..length { + state = state.wrapping_mul(1_103_515_245).wrapping_add(12345); + seq.push(bases[((state >> 16) & 3) as usize]); + } + seq +} + +fn build_genome() -> Vec { + let mut genome = lcg_seq(88888, 20_000); + for (start, end) in [INTRON, INTRON2] { + genome[start] = b'G'; + genome[start + 1] = b'T'; + genome[end - 2] = b'A'; + genome[end - 1] = b'G'; + } + genome +} + +fn write_fasta(dir: &Path, genome: &[u8]) -> PathBuf { + let path = dir.join("genome.fa"); + let mut f = fs::File::create(&path).unwrap(); + writeln!(f, ">chr1").unwrap(); + f.write_all(genome).unwrap(); + writeln!(f).unwrap(); + path +} + +fn build_index(fasta: &Path, genome_dir: &Path) { + fs::create_dir_all(genome_dir).unwrap(); + cargo_bin_cmd!("rustar-aligner") + .args([ + "--runMode", + "genomeGenerate", + "--genomeDir", + genome_dir.to_str().unwrap(), + "--genomeFastaFiles", + fasta.to_str().unwrap(), + "--genomeSAindexNbases", + "7", + ]) + .assert() + .success(); +} + +/// Reads spanning both junctions plus unspliced filler, enough of them that +/// several threads each hold some. +fn write_fastq(dir: &Path, genome: &[u8]) -> PathBuf { + let path = dir.join("reads.fq"); + let mut f = fs::File::create(&path).unwrap(); + let mut n = 0usize; + + let spliced = |f: &mut fs::File, e_end: usize, i_end: usize, overhang: usize, n: &mut usize| { + let left = &genome[e_end - overhang..e_end]; + let right = &genome[i_end..i_end + (100 - overhang)]; + let mut seq = left.to_vec(); + seq.extend_from_slice(right); + writeln!(f, "@sj{n}").unwrap(); + f.write_all(&seq).unwrap(); + writeln!(f, "\n+\n{}", "I".repeat(seq.len())).unwrap(); + *n += 1; + }; + + for overhang in 30..70 { + spliced(&mut f, INTRON.0, INTRON.1, overhang, &mut n); + spliced(&mut f, INTRON2.0, INTRON2.1, overhang, &mut n); + } + for i in 0..200usize { + let start = 8_000 + i * 40; + writeln!(f, "@u{i}").unwrap(); + f.write_all(&genome[start..start + 100]).unwrap(); + writeln!(f, "\n+\n{}", "I".repeat(100)).unwrap(); + } + let _ = (EXON1, EXON2, EXON3, EXON4); + path +} + +/// Align with `threads` threads and return the bytes of every output file that +/// carries an order. +fn align( + genome_dir: &Path, + fastq: &Path, + prefix: &str, + threads: &str, + extra: &[&str], +) -> Vec<(String, Vec)> { + let mut cmd = cargo_bin_cmd!("rustar-aligner"); + cmd.args([ + "--runMode", + "alignReads", + "--genomeDir", + genome_dir.to_str().unwrap(), + "--readFilesIn", + fastq.to_str().unwrap(), + "--runThreadN", + threads, + "--outFileNamePrefix", + prefix, + ]); + cmd.args(extra); + cmd.assert().success(); + + ["SJ.out.tab", "Aligned.out.sam"] + .iter() + .map(|name| { + let bytes = fs::read(format!("{prefix}{name}")).unwrap_or_default(); + // Drop the @PG header line: it records the command line, which + // differs by the thread count itself. + let filtered: Vec = if *name == "Aligned.out.sam" { + String::from_utf8_lossy(&bytes) + .lines() + .filter(|l| !l.starts_with("@PG")) + .collect::>() + .join("\n") + .into_bytes() + } else { + bytes + }; + ((*name).to_string(), filtered) + }) + .collect() +} + +#[test] +fn output_is_identical_at_one_and_eight_threads() { + let dir = TempDir::new().unwrap(); + let root = dir.path(); + let genome = build_genome(); + let genome_dir = root.join("genome"); + build_index(&write_fasta(root, &genome), &genome_dir); + let fq = write_fastq(root, &genome); + + let one = align( + &genome_dir, + &fq, + &format!("{}/t1_", root.display()), + "1", + &[], + ); + let eight = align( + &genome_dir, + &fq, + &format!("{}/t8_", root.display()), + "8", + &[], + ); + + // The fixture has to actually produce junctions, or this test would pass + // on two empty files. + let sj = &one[0].1; + let rows = String::from_utf8_lossy(sj).lines().count(); + println!("fixture produced {rows} SJ.out.tab rows"); + assert!( + rows >= 2, + "the fixture has to produce at least two junction rows for an order to \ + exist at all, got {rows}" + ); + + for ((name, a), (_, b)) in one.iter().zip(eight.iter()) { + assert_eq!( + String::from_utf8_lossy(a), + String::from_utf8_lossy(b), + "{name} differs between 1 and 8 threads" + ); + } +} + +#[test] +fn output_is_identical_across_two_runs_at_eight_threads() { + // Guards against the pair above agreeing only because both runs happened + // to hit the same map order. + let dir = TempDir::new().unwrap(); + let root = dir.path(); + let genome = build_genome(); + let genome_dir = root.join("genome"); + build_index(&write_fasta(root, &genome), &genome_dir); + let fq = write_fastq(root, &genome); + + let a = align( + &genome_dir, + &fq, + &format!("{}/a_", root.display()), + "8", + &[], + ); + let b = align( + &genome_dir, + &fq, + &format!("{}/b_", root.display()), + "8", + &[], + ); + + for ((name, x), (_, y)) in a.iter().zip(b.iter()) { + assert_eq!( + String::from_utf8_lossy(x), + String::from_utf8_lossy(y), + "{name} differs between two 8-thread runs" + ); + } +} + +#[test] +fn two_pass_output_is_identical_at_one_and_eight_threads() { + // Two-pass feeds pass 1's junctions back into the alignment, so any order + // escaping the junction map would show up here rather than in a single + // pass. + let dir = TempDir::new().unwrap(); + let root = dir.path(); + let genome = build_genome(); + let genome_dir = root.join("genome"); + build_index(&write_fasta(root, &genome), &genome_dir); + let fq = write_fastq(root, &genome); + + let extra = ["--twopassMode", "Basic"]; + let one = align( + &genome_dir, + &fq, + &format!("{}/tp1_", root.display()), + "1", + &extra, + ); + let eight = align( + &genome_dir, + &fq, + &format!("{}/tp8_", root.display()), + "8", + &extra, + ); + + for ((name, a), (_, b)) in one.iter().zip(eight.iter()) { + assert_eq!( + String::from_utf8_lossy(a), + String::from_utf8_lossy(b), + "two-pass {name} differs between 1 and 8 threads" + ); + } + + // The pass-1 junction file is written from the same map and is equally + // order-sensitive. + let p1 = fs::read(format!("{}/tp1__STARpass1/SJ.out.tab", root.display())).unwrap(); + let p8 = fs::read(format!("{}/tp8__STARpass1/SJ.out.tab", root.display())).unwrap(); + assert_eq!( + String::from_utf8_lossy(&p1), + String::from_utf8_lossy(&p8), + "pass-1 SJ.out.tab differs between 1 and 8 threads" + ); +}