I have tested Bismark Rust Suite v3.1.0 and compared it against Bismark v0.25.1 Perl version with methylseq pipeline. I have used the same input with the same parameters. I was hoping to have “*byte-identical*” results as reported, however, I got slightly different BAM file, which is worrying for updating production deployments and a blocker to move further to get benefit from the enhancements. The difference comes from the choice of the duplicate copy: Bismark 3.1.0 makes a real change, but it is noise, not a shift. • The aligned reads are the same, but 3.1.0 keeps a different copy from many PCR-duplicate sets. That swaps 16–44 % of each deduplicated BAM’s records. Read counts, insert sizes and overall methylation stay the same. • On the target CpGs, methylation moves by 0.08–0.10 points on average (correlation ≥ 0.99987), and by at most 11–15 points at single low-depth CpGs. • There is a small lean in all three samples: 3.1.0’s kept copies are slightly more methylated, by 0.001–0.004 points per CpG, with 0.2–1.2 % more methylated CHG/CHH calls. That is at least 20× smaller than the random part. • The lambda and pUC controls change the same way, but their gated CpG percentages don’t move. Is that related to the use of

--multicore

option? What are the potential sources of differences between the two versions? Is there a way for a “safe” upgrade? @Felix Krueger

let me look into this and get back to you

Further investigation revealed the root cause: The human

BISMARK_ALIGN

first puts each sample’s reads into the order in which Perl Bismark 0.25.1 run with

--multicore 4

aligned them, and then aligns them in that order.

deduplicate_bismark

keeps the first read of each PCR-duplicate set in alignment output order. Perl 0.25.1 with

--multicore N

gave record k of each FASTQ file to worker ((k - 1) mod N) + 1 and concatenated the workers’ BAM and unmapped-read files in worker order. 3.1.0 writes its output in input order, with

-p

and with its own

--multicore

alike, so it kept another read from many duplicate sets.

I have also just sent CLaude on the case, and this is its verdict: Exactly right, and thanks for tracking it down. The formula matches the Perl code. Reordering the input into the

--multicore 4

order is a valid way to reproduce the old dedup choices with Bowtie 2, as long as N matches what each baseline sample actually used. I wouldn’t assume it for HISAT2. For production I’d suggest re-baselining once rather than keeping the reorder step: 3.1.0's output order is the same for any core count, while Perl

--multicore N

gives a different dedup result for each N.