[jmschrei]

This optimization was found with auto-optimize, a tool developed by the Programmable Genomics Laboratory at UMass Chan Medical School: it cuts the median wall time of MarkDuplicates on a 22-million-pair DNase-seq BAM at 8 cores from 234.2 s to 25.05 s, and speeds it up 2.0x to 9.4x on 56 public BAM files, with output byte-identical to master on every input both complete. The output was checked against master in 28,294 comparisons, 28,225 of them identical, on 2,253 inputs from 50 assays (1,899 BAM, 198 SAM and 156 CRAM), among them 62 BAMs made for this check with empty BGZF blocks inserted, 14 BAMs made for this check that are cut short or have a damaged BGZF block, 10 BAMs with one BGZF defect or oddity, 14 BAMs with one record or header defect, 2 SAMs with one defect, 2 CRAMs with one defect and 428 other edge-case files from the test collection that are mostly unusual but valid; see How the output was checked.

The 22-million-pair input is ENCODE's T-cell DNase-seq BAM ENCFF377MQA (21,992,283 pairs, PE76, coordinate-sorted, one @RG added), run with -Xmx8g on 8 cores. On it, master spends 55 s reading the input and building read ends, 5 s finding duplicates, and 174 s in the second pass, where one thread decodes, marks, encodes and compresses every record. Most of the gain comes from four things: compressing BGZF blocks on several threads (items 1 and 2 take the second pass from 174 s to 41 s), decoding and validating records on several threads in both passes, writing each record from the bytes it was read from, and finding duplicates and writing the output while the first pass is still reading. Peak memory falls from 8,051 MB to 3,694 MB on that input and is lower on 54 of the 56 inputs timed (see Memory). No command-line argument is added, no default changes, and the code is plain Java 17 with no native or OS-specific calls.

It builds on #2062, which adds 54 test files for MarkDuplicates and the classes it uses and changes no source. The first five commits here are #2062's. This pull request's own change is the sixth commit: 12 files, +6,005 / -73, two of them new test classes. Once #2062 is merged, only the sixth commit remains. The bugs that those tests and the checks below found in master are reported in #2061.

Changes

In the order they were made. The time after each item is the median of three runs on the input above with every change up to that item; master takes 234.2 s.

1. Asynchronous writer (MarkDuplicates.doWork): the output writer is made with SAMFileWriterFactory.setUseAsyncIo(true), so htsjdk encodes and compresses records on its writer thread while the main thread reads and marks. Outputs that item 2's writer does not take still use it. 211.6 s.
2. Parallel BGZF compression (new ParallelBlockCompressedOutputStream and ParallelBAMFileWriter): a BGZF stream that cuts blocks exactly where BlockCompressedOutputStream cuts them (when its 65,498-byte buffer is full, and on flush()), deflates them on a pool with deflaters from the same factory at the same level, keeps the level-0 fallback for a block whose output does not fit, writes them in order and ends with the same EOF block. ParallelBAMFileWriter writes the header and records through it as BAMFileWriter does. It is used when canWrite(OUTPUT) holds: the name ends in .bam and neither an index nor an MD5 file is requested. The second pass fell from 151.8 s to 41.2 s. 101.4 s.
3. Asynchronous read in the first pass (AbstractMarkDuplicatesCommandLineProgram.openInputs(boolean, boolean)): an overload that turns on htsjdk's asynchronous BGZF read. openInputs(boolean) passes false, which keeps htsjdk's default. The first pass used it until item 12 turned it off, because htsjdk 5.0.0's asynchronous stream never returns from a read after a corrupt deflate block. No caller passes true now. 95.84 s.
4. Read-ahead thread in the first pass (new ReadAheadIterator): the input iterator runs on a background thread that hands records over in batches of 1,024, so decoding and validation overlap with building read ends. htsjdk's AsyncBufferedIterator was not used because its close() can hang: it interrupts the reading thread, htsjdk turns the interrupt into an exception and clears it, and the thread then blocks on its full queue (3 of 1,000 closes hung in a test). ReadAheadIterator.close() never interrupts. Since item 6 it is used only for inputs the parallel reader does not take. 84.58 s.
5. Paired ends out of fragSort (MarkDuplicates.buildSortedReadEndLists and generateDuplicateIndexes, new PairedReadEndKeys): the fragment pass marks only unpaired ends, and a paired end there only tells its chunk that a pair shares the chunk's key. Yet fragSort held one object per mapped primary read (44.0 million on the input above, over its 31.1 million in-memory limit, so it spilled to disk). A paired end now adds only its key (library, reference, coordinate, orientation) packed into a long, and a chunk has pairs when its key is among them. A key that cannot be packed leaves its end in fragSort, and keys that outgrow the same memory budget go to fragSort as stand-in paired ends. Peak memory fell from 8,032 to 6,289 MB. 64.76 s.
6. Parallel BAM reader (new ParallelBAMRecordIterator; an openInputs overload with a decoding-thread count): for one local BAM input on 4 or more processors, one thread inflates with htsjdk's synchronous BlockCompressedInputStream and cuts the stream into 64 KB batches of whole records, and decoding threads run BAMRecordCodec.decode and then what BAMFileReader's iterator does to each record: set the stringency, isValid with SAMUtils.processValidationErrors and the 1-based record number, and eager decoding. Records come out in file order. Like BAMFileReader, the iterator reads one record ahead, so a failure in record k is thrown by the same next() with the same exception, and LENIENT warnings print when the consumer reaches their record. The first pass fell from 20.09 to 16.17 s. 61.34 s.
7. Records written from their input bytes (new RawBlockBAMRecord; ParallelBAMFileWriter.writeFromBlock; an openInputs overload with a SAMRecordFactory): the second pass changes tags on every record (CLEAR_DT and the PG tag), so BAMRecordCodec re-encoded every field of every record. The second pass now reads RawBlockBAMRecords, which keep the variable-length block they were decoded from until the read name, cigar, sequence or qualities are set, and the writer assembles each record in one buffer: the fixed fields by direct stores, the name, cigar, sequence and qualities copied from the block, the bin computed as BAMRecordCodec computes it, and the tags through BinaryTagCodec.writeTag. Any record it cannot copy exactly (a cigar from a CG tag, out-of-range fields, a reference too long for BAI bins, a position past 2^29) goes through BAMRecordCodec. 55.39 s.
8. Batched hand-off to the writer thread (ParallelBAMFileWriter.AsyncWriter): in place of htsjdk's AbstractAsyncWriter, which costs the main thread a queue put and two atomic swaps per record, records go to the writer thread in arrays of 1,024 through a two-slot queue. Failures surface as they do from AbstractAsyncWriter. The second pass fell from 37.35 to 27.33 s. 43.31 s.
9. Parallel reader in the second pass (MarkDuplicates.doWork, ParallelBAMRecordIterator with a record factory): the second pass reads one BAM input through item 6's reader with 2 decoding threads, which make RawBlockBAMRecords. The second pass fell from 27.33 to 21.3 s. 38.19 s.
10. Duplicate-index collections sized to the input (MarkDuplicates.sortIndicesForDuplicates): htsjdk's SortingLongCollection allocates and zeroes its whole in-memory array up front, 268,435,456 longs (2 GiB) at -Xmx8g. At most one index per record of the first pass is ever added, so the size is capped at that record count. Peak memory fell from 7,715 to 5,918 MB, and the time did not change measurably. 38.39 s.
11. Read ends built on the decoding threads (MarkDuplicates.decodeRead and decodeReadEnds; an annotator in ParallelBAMRecordIterator; ReadNameParser.getReadNameRegex; an openInputs overload): for each mapped primary read, the decoding threads build the ReadEndsForMarkDuplicates that buildReadEnds would (unclipped 5' coordinate, orientation, score, mate reference, and tile, x and y from a ReadNameParser of their own), build the pair map's key (read group ID and read name, with its hash computed), and read the PG tag. The main thread sets the library id and file index in file order. If any lookup throws, the main thread repeats all of them with master's code, so the same exception surfaces at the same record. The first pass uses min(4, processors - 2) decoding threads. It fell from 15.5 to 10.04 s. 33.36 s.
12. Damaged BGZF files read as BAMFileReader reads them (ParallelBAMRecordIterator, MarkDuplicates.buildSortedReadEndLists): an empty BGZF block, or one whose ISIZE is 0, ends the records only where it ends BAMFileReader's (cleanly at a record's start, with "Premature EOF" where BAMRecordCodec.decode starts one of its other reads) and is read past elsewhere. An IOException from the seek to the first record is wrapped as BAMFileReader wraps it, and the first pass no longer uses htsjdk's asynchronous BGZF read (item 3). On 78 BAMs with empty, zero-length or corrupt blocks, the reader matches htsjdk's on every one, and MarkDuplicates exits with master's code and first exception line in 16 of 16 cases. 33.09 s.
13. One validation per record (ParallelBAMRecordIterator.open and an openInputs overload with a count of validated records; MarkDuplicates.doWork): the first pass has already validated each record it read, at the same stringency, so the second pass's reader skips isValid for those and decodes the tags and then the cigar at SILENT, which is what validation would have decoded. Later records are validated as before. See Behaviour changes for the LENIENT warnings. 31.95 s.
14. Fixed fields decoded from the bytes (ParallelBAMRecordIterator.createRecord): in a batch that is not the last, whose records the splitter has already checked, records are made from little-endian loads with the same conversions and factory call as BAMRecordCodec.decode, rather than through 13 stream reads each. The last batch, which can end in a partial record or a stored failure, still goes through BAMRecordCodec. 31.13 s.
15. Tags left undecoded in the second pass (RawBlockBAMRecord.deferTags and writeTags; ParallelBAMRecordIterator.decodeValidated; ParallelBAMFileWriter.writeFromBlock): a validated record whose tags and cigar cannot fail to decode (every cigar operator known, the tags filling the block exactly, no CG tag, no repeated tag, no H value, no B subtype outside c, C, s, S, i, I and f) keeps them as bytes. getAttribute reads from the bytes, setAttribute with null or a String is recorded as an edit, and anything else decodes everything at SILENT and replays the edits. The writer merges the original tags with the edits as BinaryTagCodec would write them. 30.57 s.
16. Deflater threads (ParallelBAMFileWriter.defaultThreads): max(1, min(p - p/4, 8)) for p available processors, 6 on 8 processors, leaving a quarter to the threads that read, decode and encode. With 4 or 5 deflaters on 8 processors, deflation bounds the second pass. 29.87 s.
17. Completed pairs as rows (new PairedReadEndRows; MarkDuplicates.buildSortedReadEndLists and generateDuplicateIndexes): each completed pair is a row of 7 longs in a 64-byte slot, in chunks of up to 2^20 rows, rather than an object in pairSort, so young collections no longer copy 22 million objects. sort() orders the rows as ReadEndsMDComparator orders the ends (packed prefixes sorted with Arrays.parallelSort, ties ordered by the full comparison) and marks where duplicate sets start, and a set of one pair is only counted. An end the rows cannot hold exactly (an int the comparator would subtract outside [-2^30, 2^30), or more pairs than the in-memory limit) moves every row into master's pairSort. G1 pauses in the first pass fell from 642 ms to 122 ms in one logged run, and peak memory from 5,643 to 5,147 MB. 29.18 s.
18. Pair pass during the first pass (MarkDuplicates.PairPassDuringRead, markDuplicatePairs): for one coordinate-sorted input, the main thread also hands completed pairs to a helper thread in batches of 4,096, each with a frontier 1,000 bp below the alignment start of the last read. Once a duplicate set's read2 end lies below a frontier, every pair of the set has been read, and the helper handles the set as handleChunk does, into its own indices and histograms, which are added after the first pass. A read below the frontier, a pair the rows refuse, too many waiting pairs, or a failure abandons the helper's result, and the pair pass runs on the rows as in item 17. The time between the passes fell from 1.13 to 0.10 s, and peak memory from 5,147 to 4,345 MB. 28.36 s.
19. The writing thread helps deflate (ParallelBlockCompressedOutputStream.writeOldest and deflateQueuedBlock): while it waits for the oldest block, the thread that writes runs deflate tasks no worker has started. Their futures are already queued in order, so blocks are still written in order. 27.99 s.
20. Keys sorted on a helper thread (PairedReadEndKeys): keys are added into chunks of 2^20 that a sorter thread sorts, de-duplicates and places into ordered runs, so doneAdding sorts only the last chunk.
21. One read-group lookup per read (MarkDuplicates.decodeRead and ReadGroupSlots): the decoding threads look up each mapped primary read's read group once, and the main thread takes its library id from an array indexed by the read group's place in the header.
22. Single-threaded progress logger (MarkDuplicates.SingleThreadProgressLogger): logs the same lines at the same records as htsjdk's ProgressLogger, without its lock, and reads the read name only for a line that prints it. Items 20 to 22 were measured together: peak memory fell from 4,347 to 3,764 MB. 27.18 s.
23. Second pass during the first pass (MarkDuplicates.SecondPassDuringRead, startSecondPass and writtenDuringRead; ParallelBAMRecordIterator.ValidatedRecords; a create-new option and a blocks-in-flight setter in ParallelBAMFileWriter; ParallelBlockCompressedOutputStream.setBlocksInFlightLimit): the second pass runs on a thread of its own while the first pass reads, and writes record i only once every duplicate flag below record i+1 is decided. Item 18's helper now decides every flag: those of pairs as before, and those of unpaired ends a chunk at a time, once the chunk's 5' position lies below a frontier, as the fragment pass would; the main thread hands it the unpaired ends and the keys of paired ends for this. A read still waiting for its mate when the frontier passes the mate's alignment start is taken to have no mate, and so not to be a duplicate. The second pass writes a new hidden file beside OUTPUT and moves it onto OUTPUT at the end. It holds back its log lines and releases them after the first pass, in master's order. While the first pass reads, 96 BGZF blocks are kept in flight instead of 24. When any condition in Behaviour changes fails, the thread stops, its file is deleted, and the second pass runs after the first as in master. 25.70 s.
24. Fallback memory (PairedReadEndKeys.splitRuns and its overflow path; MarkDuplicates.generateDuplicateIndexes; new PairedReadEndKeysTest and PairedReadEndRowsTest): before master's order allocates the duplicate indices, the sorted key arrays are split into arrays of at most 32,768 keys (256 KiB), and the rows are released once the helper has handled every pair. The key overflow path splits the same way before it hands keys over. See Memory. The input above does not take these paths. 25.05 s.

Tests and CI

The tests from #2062 and the two new test classes pass on this branch. The repository's Build and Test workflow was run locally, step by step, on this commit with OpenJDK 17:

Job	Result
./gradlew compileJava, shadowJar	passed
legacyTest and barclayTest filtered to picard.sam.markduplicates.*, picard.sam.util.*, picard.sam.DuplicationMetrics* and picard.sam.FlowBasedDuplicationMetricsTest	2,712 passed in each
legacyTest (which jacocoTestReport runs), every test	7,140 passed
barclayTest, every test	7,139 passed
./gradlew currentJar, then java -jar build/libs/picard.jar MarkDuplicates -I testdata/picard/sam/aligned_queryname_sorted.bam -O out.bam --METRICS_FILE out.metrics	passed; with the @PG version pinned, out.bam and the metrics are identical to master's

Every test passed on this commit, as on master (5,466 and 5,465 tests there). In earlier runs launched another way on the same machine, five tests that pipe data through standard input or output (FastqToSamTest.testStreamInput, PipedDataTest.testPipedInput[0..2] and UpdateVcfSequenceDictionaryTest.testCaptureStdout) failed on master and this branch alike; launched as above, they pass on both. The Cloud Tests workflow needs Google Cloud credentials and was not run, and neither was the Trivy scan.

Behaviour changes

There are four.

• Partial output after a failure. When MarkDuplicates fails while writing, both builds exit 1 with the same exception, but the partial file left behind differs. In the comparisons below this happened on inputs whose records are not in the order their header states ("Alignments added out of order") and once on a BAM index failure. master writes on the main thread and leaves a flushed partial file. This branch writes on a writer thread; its failure surfaces on a later addAlignment, after which close() does not close the underlying writer, as in htsjdk's AbstractAsyncWriter. The partial file is then often empty, or a shorter SAM, and which one is left varies from run to run: on one such input with REMOVE_SEQUENCING_DUPLICATES=true TAGGING_POLICY=All, this branch left an empty file in 5 of 8 runs and master's 116,043-byte file in 3, while master left that file in 4 of 4. 69 of the 28,294 comparisons below differ in this way and in no other.
• LENIENT validation warnings. master validates every record in both passes, so at VALIDATION_STRINGENCY=LENIENT it prints each warning twice for an invalid record that both passes read. Where the second pass uses the parallel reader, this branch does not validate those records again (item 13) and prints the warning once. In the comparisons below, the number of warnings differed in 373 jobs: 360 printed fewer (245 of them exactly half as many) and 13 printed more, 12 of them on three SAM inputs that both builds then reject with the same exception (up to 62 warnings against 2) and one on a BAM (3 against 2). In a separate check, a SAM input and a run with two inputs printed one record's warning 4 times against master's 2. We did not establish why these cases print more. Exit codes, exceptions and outputs are the same in all of them.
• When the passes overlap. The second pass runs during the first (item 23) only when all of these hold: MarkDuplicates itself rather than a subclass, without barcode tags or flow mode; one INPUT, a local coordinate-sorted BAM file; TAG_DUPLICATE_SET_MEMBERS=false; READ_NAME_REGEX null or the standard optical duplicate finder; at least 4 available processors and a heap of at least 1 GiB; and an OUTPUT that does not exist yet, whose name ends in .bam, without CREATE_INDEX or CREATE_MD5_FILE. It is also abandoned during the run on a read or unpaired end below the frontier, a pair or key that cannot be held, too many pairs, ends, keys or reads waiting, a read with a mapped mate whose RNEXT is *, a read taken to have no mate that then finds it, a PG ID the output header lacks, or a failure in the first pass. Otherwise the passes run in master's order with every other change; on the input above with an existing OUTPUT, that took 33.32 to 34.47 s in three runs against master's 238.6 to 243.3 s. While the passes overlap, the output is written to a hidden file .<OUTPUT name>.<random hex>.tmp in OUTPUT's directory and moved onto OUTPUT at the end (an atomic move where the file system allows one), so OUTPUT appears only when the run finishes, and a process that is killed can leave the hidden file behind. No command-line argument chooses the path, because Picard writes every argument into the output's @PG CL, so a new one would change the output bytes.
• Logs and stack traces. The "Will retain up to N duplicate indices" line reports the capped size (item 10), the order of "Large duplicate set" lines can change (item 18), and an exception may carry a stack trace from a reading or writing thread rather than the main thread, with the same class and message.

Usage benchmarks

Machine: AMD EPYC 9575F, 2 sockets x 64 cores (256 hardware threads, 2 NUMA nodes), 1 TB RAM, Ubuntu 24.04.5, Linux 7.0, OpenJDK 17.0.18, htsjdk 5.0.0 with Picard's default Intel deflater and inflater. Both builds are ./gradlew shadowJar jars: master (with #2062's test files, which change no source) and this branch.

Inputs: no synthetic data. Every input is public sequencing data from a test collection: one large subsample per assay (46) and 10 full-depth experiments, two of them limited to chr22. Paired-end files without a read group got one @RG with samtools addreplacerg, because master throws a NullPointerException on paired reads without one (#2061). One whole-genome subsample is left out of the table for its license.

How it was timed:

• Isolation: one configuration at a time, pinned with taskset to the 8 physical cores of one L3 group (CPUs 24-31), to which no other process was pinned. Other jobs ran on the rest of the machine; the load average at the run starts was 61.5 to 740.8 (median 151.3).
• Repeats: three runs per build, alternating which build went first. The table gives the median.
• What the time includes: wall time from process start to exit, JVM start-up and closing the output included, with -Xmx8g and VALIDATION_STRINGENCY=LENIENT, so that a file with a few invalid records runs to the end; every record is still validated. The input and the jar were in the page cache, and the output went to local disk.
• Peak memory: the process's maximum resident set (/usr/bin/time %M).
• Output check: whether out.bam was byte-identical and the metrics identical apart from their timestamp and command-line lines, with --PROGRAM_GROUP_VERSION set to the same value for both builds so that the @PG line matches.

Input class	Input	Threads	master (s)	This branch (s)	Speedup	master peak (MB)	This branch peak (MB)	Output
benchmark input	DNase-seq T cell (ENCODE ENCFF377MQA), PE76, 21,992,283 pairs, @rg added	8	240.38	25.48	9.44x	7,813	3,674	identical
full depth	ChIP-seq lung input control, PE76, 29,922,312 alignments, @rg added	8	185.30	21.05	8.80x	7,448	2,955	identical
full depth	csRNA-seq K562, SE, 28,950,388 alignment records	8	72.08	8.60	8.38x	4,965	3,229	identical
full depth	CUT&Tag K562 H3K27me3, PE25, 4,120,414 pairs	8	39.63	5.51	7.19x	5,708	1,256	identical
full depth	eCLIP K562 RBFOX2, PE, 20,298,696 alignment records	8	89.14	10.03	8.89x	6,347	2,029	identical
full depth	MNase-seq K562, PE50, 22,767,173 pairs, @rg added	8	267.00	29.89	8.93x	8,133	3,831	identical
full depth	mNET-seq K562, PE, chr22 only, 4,691,806 pairs, @rg added	8	39.23	5.65	6.94x	5,720	2,268	identical
full depth	PRO-seq K562, SE, 191,353,062 reads	8	1,021.28	152.53	6.70x	4,462	6,203	identical
full depth	SLAM-seq K562, SE, 42,031,760 reads	8	198.70	31.42	6.32x	3,748	6,299	identical
full depth	WGS tumour HCC1395, PE150, chr22 only, 19,906,411 alignments	8	121.54	16.29	7.46x	6,298	2,493	identical
large subsample	10x Multiome gene expression, PBMC, 434,740 reads	8	4.57	1.27	3.59x	3,676	1,127	identical
large subsample	10x scATAC-seq PBMC, PE50, 981,668 alignment records	8	6.83	1.67	4.09x	4,117	1,140	identical
large subsample	10x scRNA-seq PBMC, 3' v3, 862,498 reads	8	8.18	1.58	5.19x	4,913	1,264	identical
large subsample	Amplicon SMF mESC, PE, 374,436 alignments, @rg added	8	8.08	1.97	4.09x	3,674	1,142	identical
large subsample	Ancient DNA ORT16, SE collapsed, 249,957 alignments	8	2.62	0.92	2.84x	3,175	822	identical
large subsample	ATAC-seq A549, SE51, 400,286 reads	8	3.37	0.97	3.48x	3,594	765	identical
large subsample	BG4 CUT&Tag K562, PE36, 500,174 pairs, @rg added	8	6.28	1.58	3.98x	4,093	874	identical
large subsample	CAGE GM12878 nucleolus, SE36, 1,004,318 alignments	8	5.12	1.22	4.19x	3,659	1,109	identical
large subsample	Capture-C mESC, PE50, 316,648 alignments, @rg added	8	3.37	1.32	2.56x	3,459	956	identical
large subsample	cfDNA plasma, PE39, 474,812 alignments	8	4.17	1.27	3.28x	3,605	1,102	identical
large subsample	ChIA-PET K562 CTCF, PE, 499,845 alignments, @rg added	8	7.38	2.37	3.11x	3,669	1,066	identical
large subsample	ChIP-exo yeast Reb1, PE40/36, 499,958 pairs, @rg added	8	6.42	1.53	4.21x	4,362	839	identical
large subsample	ChIP-nexus mESC Nanog, SE76, 1,001,291 reads	8	6.98	1.38	5.07x	4,388	1,230	identical
large subsample	ChIP-seq HepG2 input control, SE100, 399,245 alignments	8	3.97	0.97	4.08x	3,623	794	identical
large subsample	csRNA-seq K562, SE, 500,852 reads	8	3.37	1.02	3.31x	3,600	1,104	identical
large subsample	CUT&RUN K562 IgG, PE51, 499,779 pairs	8	6.93	1.52	4.55x	4,261	1,112	identical
large subsample	CUT&Tag K562 H3K27me3, PE25, 500,128 pairs	8	6.63	1.47	4.51x	4,586	1,090	identical
large subsample	DMS-MaPseq HEK293T, SE48, 500,228 alignments	8	2.92	0.92	3.17x	3,531	959	identical
large subsample	DNase-seq heart, SE36, 401,154 reads	8	3.42	0.97	3.54x	3,582	678	identical
large subsample	DRIP-seq NT2 RNase H control, SE100, 499,369 reads	8	5.72	1.17	4.88x	3,642	1,080	identical
large subsample	eCLIP K562 RBFOX2, PE, 505,988 alignment records	8	3.87	1.12	3.45x	3,628	712	identical
large subsample	END-seq MCF-7, SE75, 500,357 reads	8	4.57	1.17	3.91x	3,630	832	identical
large subsample	Exome HG002, PE151, 500,515 alignments	8	5.83	1.37	4.24x	3,669	1,094	identical
large subsample	Fiber-seq K562, long reads, 99,685 alignments	8	76.64	11.02	6.95x	3,772	2,259	identical
large subsample	GLORI HEK293T, SE, 500,793 alignments	8	4.12	1.07	3.85x	3,602	1,080	identical
large subsample	HiChIP K562 H3K27ac, PE76, 372,492 alignments, @rg added	8	4.53	1.57	2.88x	3,578	932	identical
large subsample	iCLIP HEK293 PTBP1, SE, 500,260 alignments	8	4.28	1.12	3.82x	3,615	874	identical
large subsample	MNase-seq K562, PE50, 1,000,601 pairs, @rg added	8	13.80	2.33	5.92x	4,552	880	identical
large subsample	mNET-seq K562, PE, 500,195 pairs, @rg added	8	6.42	1.52	4.22x	4,147	849	identical
large subsample	NOMe-seq K562, PE, 500,126 alignments, @rg added	8	7.28	1.43	5.10x	3,617	1,091	identical
large subsample	PacBio HiFi HG002, long reads, 104,103 alignments	8	83.65	17.68	4.73x	3,659	1,693	identical
large subsample	PolyA-seq GM12878, SE, 515,751 alignments	8	4.92	1.22	4.03x	3,620	1,063	identical
large subsample	Pore-C GM12878, nanopore, 217,901 alignments	8	5.77	1.28	4.52x	3,641	1,078	identical
large subsample	PRO-cap K562, PE, 568,462 pairs, @rg added	8	7.13	1.68	4.26x	4,928	1,202	identical
large subsample	PRO-seq K562, SE, 500,206 reads	8	4.82	1.07	4.50x	3,634	1,072	identical
large subsample	Repli-seq K562 late S, SE, 577,049 reads	8	4.32	1.12	3.86x	3,809	1,022	identical
large subsample	Ribo-seq K562, SE51, 958,882 alignment records	8	4.22	1.17	3.61x	3,687	1,139	identical
large subsample	RNA-seq K562 poly(A), SE76, 585,411 alignment records	8	4.73	1.27	3.72x	3,637	946	identical
large subsample	SARS-CoV-2 ARTIC amplicons, long reads, 9,382 alignments	8	1.22	0.61	1.98x	2,642	383	identical
large subsample	scBS-seq mESC, SE, 399,766 alignments	8	6.32	1.27	4.98x	3,628	1,077	identical
large subsample	SLAM-seq K562, SE, 499,144 reads	8	4.32	1.22	3.55x	3,674	1,137	identical
large subsample	STARR-seq K562 DNA input, PE150, 398,913 pairs, @rg added	8	10.48	1.77	5.92x	4,096	1,132	identical
large subsample	Start-seq T-47D, PE, 500,291 pairs, @rg added	8	5.88	1.37	4.29x	4,457	1,117	identical
large subsample	TT-seq K562, PE76, 500,287 pairs, @rg added	8	8.58	1.68	5.11x	4,127	1,107	identical
large subsample	Visium breast cancer, SE50, 576,233 reads	8	4.77	1.68	2.85x	3,994	1,068	identical
large subsample	WGS tumour HCC1395, PE150, 513,272 alignments	8	4.82	1.32	3.65x	3,686	1,110	identical

The output was byte-identical in every configuration. On the 22-million-pair input, a fresh clone of master with this change merged, both built with ./gradlew shadowJar, gave 238.6 to 243.3 s for master and 25.11 to 25.95 s for this branch in three runs each, and a paired comparison of master with an earlier version of this change (before item 24) found the change faster in 9 of 9 pairs (median difference -213.4 s; exact two-sided sign test, p = 0.0039). Over the 56 rows, the speedup has a geometric mean of 4.49x and a median of 4.21x, and the summed time falls from 2,678 s to 395 s. Over the 10 full-depth rows the geometric mean is 7.84x (6.32x to 9.44x), and over the 46 subsamples 3.98x (1.98x to 6.95x).

How the output was checked

The usage benchmarks above check speed and memory. Separately, we ran master and this branch side by side on every alignment file in the same test collection:

• Scale: 28,294 comparisons over 2,253 inputs. The inputs are the 1,579 BAM, SAM and CRAM files of the collection (50 assays, tiny to full depth, 150 of them full-depth, 523.6 GB in all, 456 of them edge cases or malformed), 519 copies with @RG added, 49 queryname-sorted copies, 30 copies with MC tags added by samtools fixmate -m, and 78 BAMs with empty, zero-length or corrupt BGZF blocks made for this check. Two @RG copies could not be made, so 26 jobs were not compared.
• What was compared: the exit code; after a success every file the run left (out.bam or out.sam, .bai, .md5) by SHA-256, and the metrics without their timestamp and command-line lines; after a failure the first exception (class and message) and the partial output file. Every path was relative and --PROGRAM_GROUP_VERSION was the same for both builds, so @PG lines are compared too.
• Builds: the same two jars as above, with -Xmx8g (-Xmx32g for both where master ran out of memory at 8g) on 8 cores, except where noted. master's side of each comparison was recorded in an earlier run of these comparisons and run again wherever the two builds differed.
• Invalid input: a job whose master run stops at validation was run again with VALIDATION_STRINGENCY=LENIENT.

What ran

Tool and options	Inputs	Comparisons	Identical	Identical failure	Differ
MarkDuplicates, default options	every input	2,359	1,614	739	6
MarkDuplicates with each of REMOVE_DUPLICATES; REMOVE_SEQUENCING_DUPLICATES with TAGGING_POLICY=All; DUPLICATE_SCORING_STRATEGY TOTAL_MAPPED_REFERENCE_LENGTH and RANDOM; READ_NAME_REGEX=null with OPTICAL_DUPLICATE_PIXEL_DISTANCE=2500; forced spilling (SORTING_COLLECTION_SIZE_RATIO=0.00001, MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=2); CREATE_INDEX with CREATE_MD5_FILE; an existing output file; .sam output; TAG_DUPLICATE_SET_MEMBERS (queryname inputs); BARCODE_TAG=CB with MOLECULAR_IDENTIFIER_TAG=MI (single-cell inputs); two INPUTs (16 pairs of files)	every input but full depth	19,231	13,045	6,123	63
MarkDuplicatesWithMateCigar, default options	as above, plus the MC-tagged copies	2,143	1,604	539	0
SimpleMarkDuplicatesWithMateCigar, default options	as above	2,143	1,471	672	0
UmiAwareMarkDuplicatesWithMateCigar, UMI_METRICS_FILE added	as above	2,143	5	2,138	0
MarkDuplicates at LENIENT on 1, 2, 4, 8 and 16 cores, against master on 1 core	one large BAM per assay and 8 full-depth BAMs	275	275	0	0
all		28,294	18,014	10,211	69

"Identical failure" means both builds failed with the same exit code, the same first exception and the same partial output. The commonest is the NullPointerException that master throws on paired reads without a read group (4,746 comparisons).

Causes of every difference

Cause	Comparisons	Evidence
The partial output left after the writer rejects records with "Alignments added out of order", on 10 inputs: three name-sorted copies, three PacBio RNA files and four edge-case files	68	same exit code and exception in both builds; the partial file varies between runs of this branch (Behaviour changes), and in 11 of these jobs an earlier version of this change left master's partial file
The partial output after "Exception when processing alignment for BAM index" (CREATE_INDEX at LENIENT, hg19 records under an hg38 header)	1	same exit code and exception in both builds

Apart from these, every exit code, exception and output file is the same.

Controls

• master against itself: 334 comparisons on 92 inputs, all identical.
• A planted fault: one > changed to >= in markDuplicatePairs, so that the last of two tied pairs is kept: 23 of 194 comparisons differ, on 10 inputs.

Threads

The new code uses more threads, and Runtime.availableProcessors() sets how many, so taskset and container CPU limits are honoured:

• Reading: one inflating thread and min(4, p - 2) decoding threads in the first pass and 2 in the second, from 4 available processors up. Below 4, the first pass reads on one read-ahead thread and the second on the main thread.
• Writing: one writer thread and max(1, min(p - p/4, 8)) deflating threads.
• Duplicate finding: a helper for the pair pass, a sorter for the paired-end keys and a thread for the overlapping second pass, under the conditions in Behaviour changes.

The output does not depend on the thread count: the 275 comparisons at 1 to 16 cores above are all identical. CPU time falls with the wall time: on the input above, from 250.8 to 184.3 CPU-seconds (median of three runs at 8 cores). On 1 core, a paired comparison on a 400,333-pair subsample of the same input found an earlier version of this change faster in 9 of 9 pairs (median difference -0.657 s on about 7.9 s; exact two-sided sign test, p = 0.0039).

Memory

Peak resident memory is lower on 54 of the 56 inputs timed, with a median of 0.29 times master's. On the 22-million-pair input, the steps that lowered it were items 5, 10, 17, 18 and 20 to 22 (numbers in Changes). It is higher on two single-end full-depth inputs: SLAM-seq (42,031,760 reads), 3,748 MB against 6,299 MB (1.68x), and PRO-seq (191,353,062 reads), 4,462 MB against 6,203 MB (1.39x). We did not establish which change raises it there.

Above about maxMemory x SORTING_COLLECTION_SIZE_RATIO (0.25 by default) / 69 read pairs (31 million at -Xmx8g), the rows refuse a pair and the code falls back to master's order (item 17). An earlier version of this change ran out of memory there once: on a STARR-seq DNA input (ENCODE ENCFF807BAQ with @RG added, 113.6 million pairs, 18.2 GB) at -Xmx8g, it threw OutOfMemoryError in sortIndicesForDuplicates, where master finished. -Xlog:gc+humongous showed why. On the fallback, the duplicate-index SortingLongCollection needs one contiguous 1.8 GB long[] (434 regions of 4 MiB), while the sorted key arrays of PairedReadEndKeys, about 8 MiB each, were still live: G1 allocates such arrays as humongous objects and never moves them, and 210 of them lay spread through the heap. Item 24 copies the keys into arrays of at most 32,768 keys (256 KiB, under half of G1's smallest region) and releases the rows before that allocation; the full collection after it left 465 humongous regions against 884 before the fix. The failure did not recur in two later runs of the earlier version (at -Xmx8g and -Xmx7g), so the fix rests on the heap layout above and on a check that forces each fallback path at a reduced heap and compares the output with master's at the same heap:

Input and heap	Path	This branch	master	Output
T-cell DNase-seq, chr1 to chr3, 1 GB	the rows refuse a pair after 3,890,368	12.8 s	55.2 s	identical
T-cell DNase-seq above, 1 GB	the rows refuse; 30,103,688 keys handed to fragSort	46.6 s	237.4 s	identical
ChIA-PET ENCFF860BVR, chr10 to chr12, @RG added, 1 GB	524,288 reads wait for their mates; the second pass stops, the pair pass finishes, the rows are released	6.2 s	35.9 s	identical
STARR-seq above, 7 GB	the rows refuse after 27,232,582	345.4 s, peak 7,585 MB	1,904.1 s	identical
STARR-seq above, 8 GB	the rows refuse after 31,122,951	361.2 s, peak 8,378 MB	1,896.9 s	identical

"Identical" covers out.bam, the metrics without their header lines, and the files left behind. Large inputs that take the fallback still run about 5x faster than master.

What this testing does not cover

• Platforms: only x86-64 Linux (Ubuntu 24.04), one AMD Zen 5 CPU and OpenJDK 17.0.18. We tested no macOS, Windows, ARM machine or other JDK; Picard's CI runs on Ubuntu with Temurin 17. The code has no OS-specific calls.
• Timing: only at 8 cores and -Xmx8g, apart from the 1-core comparison and the reduced heaps above. The output was compared at 1, 2, 4, 8 and 16 cores, but no other core count was timed, so item 16's deflater count is untimed on 4 to 7 and 9 or more processors.
• Options not in the comparisons: flow mode, COMPRESSION_LEVEL and USE_JDK_DEFLATER. Flow mode is covered by [auto-optimize]: Add unit tests for MarkDuplicates and related tools #2062's tests, and the compression options (levels 0, 1 and 9, and the JDK deflater) were compared with master on one small input only. With the Intel deflater at COMPRESSION_LEVEL=1, master itself is not reproducible from run to run, so neither build can be held byte-identical there.
• The paired comparisons ran before item 24, which changes only the fallback paths.

Bugs in master found along the way

The bugs that #2062's tests and these checks found in master are in #2061, and those that live in htsjdk will be reported to samtools/htsjdk and linked here. This pull request fixes none of them. Among them:

• Paired reads without a read group throw a NullPointerException in MarkDuplicates (4,746 of the comparisons above fail this way in both builds).
• MAX_OPTICAL_DUPLICATE_SET_SIZE=-1, which the documentation says disables the check, turns optical duplicate detection off.
• With READ_NAME_REGEX=null, the duplicate set size histograms miss every set of two or more pairs.
• In htsjdk 5.0.0, which Picard uses, an empty BGZF block exactly between two records reads as the end of the file, so MarkDuplicates exits 0 having written 2,797 of 5,594 records; both builds write the same file. Read past an empty BGZF block that is not the last block samtools/htsjdk#1813 fixed this in htsjdk 6.0.0, so it is not filed again. MarkDuplicates does not check BGZF CRC32s, which the htsjdk report will cover.
• htsjdk 5.0.0's AsyncBlockCompressedInputStream never returns from a read after a corrupt deflate block (item 3).

🤖 Generated with Claude Code

[jmschrei]

Hello! The above message was written by Claude, but this comment is written by me (though based on my other auto-optimization PRs).

I have been developing a tool, called [auto-optimize], which basically applies the auto-research loop to code to improve speed and memory usage of code, but with extensive robustness checks to try to ensure nothing subtly breaks. I have been applying [auto-optimize] across the genomics ecosystem, starting with my own tools but now extending to other important ecosystem tools, with the hope of making things broadly better. Most recently, I've applied it to MACS3, Bowtie2, chromap, and bedtools.

In an effort to be more robust, the process first augments the unit tests to make sure that nothing breaks during the search. See #2062 (feel free to merge that independently of any choice you make here). Any iteration in our loop that broke any unit test was automatically discarded, and we used a set of external validations that were not directly optimized against each iteration to make sure we got identical results. At the end, we ran a bunch of example cases across assays and data type, covering every command-line flag, combinations, etc. Importantly, we also created a set of malformed data to make sure that not only did you get identical outputs on good data, but that it failed the same way on bad data so that we didn't just get rid of safeties that were necessary but slow.

If you look at the results below, you'll see that there were several big improvements, but also a collection of small improvements that all collectively led to the improved speed.

In an effort to be consistent with https://rewrites.bio/, I'd like this to be merged into this repository so the original developers and maintainers can continue to get credit for their work (though I don't think that's really in question...). However, as a maintainer myself, I understand being skeptical of merging a PR of this size, particularly when the submitter can't explain everything in it (a consequence of so many changes being made at a low level). There are several paths forward: (1) you can merge the PR as-is, (2) you can push your own commits on top of it to revert any changes you dislike or built upon them, (3) you can cherry-pick out any individual or set of improvements you like and merge them as separate PRs (as Ben is doing with bowtie2 right now), (4) you can ignore this entirely. We're happy to share any logs or information that would help you!

[serge2016]

Dear @jmschrei, could you also try to speedup GATK MarkDuplicatesSpark (https://github.com/broadinstitute/gatk/blob/master/src/main/java/org/broadinstitute/hellbender/tools/spark/transforms/markduplicates/MarkDuplicatesSpark.java)?

[tfenne]

Thanks for this work @jmschrei. I just wanted to respond here and say that this is both impressive, and unfortunately unlikely to go anywhere in the near future due to i) the size of the PR, and ii) the lack of any real resourcing to maintain Picard.

That said, this is quite impressive, and I'm tempted to mine this PR for certain things that could go upstream into HTSJDK which Picard relies upon and I've been trying my best to resuscitate development on lately. E.g. a lot of the parallelization work for reading/writing BAM files could probably flow upstream and be made more broadly available to both Picard tools and other tools that use HTSJDK. Is that something you would welcome and/or want to participate in?

[jmschrei]

Hi @tfenne. Totally understand. Regardless, this is here for anyone to use if they'd like and for developers to pick apart. Feel free to take any of the optimizations found here and merge/use them wherever you'd like. Credit pointing to this PR would be appreciated if you do so but tbh it's not required. Alternatively, if you have specific packages or tooling you think would welcome a PR of this type (even if just to cherry-pick out individual ideas), I'm happy to run auto-optimize on those.