CRAM 3.1 write support: full codecs, profiles, optimisation, and tests

Builds out CRAM 3.1 write support on top of cnorman's initial naive
profile, taking htsjdk to feature parity with samtools/htslib for
writing CRAM 3.1 and matching samtools on cross-implementation fidelity
across a broad range of real-world data.

CRAM 3.1 codecs and supporting infrastructure
---------------------------------------------
- Full CRAM 3.1 codec set: rANS Nx16 (Order 0/1, RLE, Stripe, Pack,
  Cat), FQZComp (quality scores), Range adaptive arithmetic coder, and
  the Name Tokeniser.
- New compression profile system (FAST, NORMAL, SMALL, ARCHIVE),
  matching htslib's `--output-fmt-option fast/small/archive`. Each
  profile picks per-DataSeries compressors and (for SMALL/ARCHIVE)
  trial-compression candidate sets.
- TrialCompressor: a wrapper that tries multiple compressors per block
  and keeps the smallest output. Replaces the ad-hoc tag triple-
  compression path with a uniform, profile-driven mechanism.
- DataSeries content IDs aligned with htslib so CRAM dumps from the
  two implementations are directly comparable.

Codec performance work
----------------------
- rANS codecs reworked to use a `byte[]` API and backwards-write,
  removing per-byte stream overhead and several layers of indirection.
- GzipCodec uses Deflater/Inflater directly instead of going through
  GZIPInputStream/OutputStream, avoiding gzip framing overhead per
  block.
- Name Tokeniser encoder: hand-written tokenisation replaces regex,
  per-type flags + STRIPE + duplicate-stream elimination + all-MATCH
  fast path significantly improve both speed and ratio.
- FQZComp / Range coder / rANS encoder hot paths tightened.
- Pooled RANSNx16Decode instance in the Name Tokeniser to avoid
  reallocating per call.
- Replaced ByteArrayInputStream/OutputStream with unsynchronized
  CRAMByteReader/Writer to remove monitor overhead in tight loops.
- Cached SAMTag key metadata to eliminate per-record String allocation
  during CRAM decode.
- Fused read-base restoration, CIGAR building, and NM/MD computation
  into a single pass instead of three.
- jlibdeflate compatibility fix on the ByteBuffer path.
- Roughly 15% faster encoding overall vs the pre-optimisation state,
  with read decoding gains in line.

Correctness fixes found during cross-implementation validation
--------------------------------------------------------------
- TLEN is now computed using htslib's coordinate-extent rule
  (max(end1,end2) - min(start1,start2) + 1, signed by leftmost
  position) for CRAM-specific compute, rather than the SAM 5'-to-5'
  rule. Without this, every read with a supplementary alignment
  mismatched on TLEN through any CRAM round-trip.
- CompressionHeader serialisation now uses a growable
  ByteArrayOutputStream rather than a fixed 100 KB ByteBuffer for the
  preservation map and tag encoding map. The TD tag dictionary alone
  can exceed 100 KB for rich tag sets (PacBio/Ultima flow-space, ONT
  modified-base tags) under the archive profile, where it would
  previously throw BufferOverflowException.
- Slice flush now uses a dual record/bases threshold matching htslib's
  `seqs_per_slice` AND `bases_per_slice` (default = readsPerSlice *
  500). Without this, archive-profile encoding of long-read data
  (PacBio HiFi 15 kb, ONT 30 kb+) would buffer ~1+ GB of quality
  scores per slice and OOM the FQZComp encoder.
- Strip NM/MD on encode and regenerate them on decode (matching
  htslib's default `store_nm=0`/`store_md=0` behaviour). Implemented
  attached mate pairs to align the in-memory representation with the
  spec.
- Restored CIGAR reconstruction when SEQ is `*` (CF_UNKNOWN_BASES).
- Fixed crash when reading containers that contain no slices.
- Removed unnecessary content digest tags from CRAM slice headers
  (htslib doesn't write them either; htsjdk's verification on read
  was overly strict).

Tools
-----
- `CramConverter`: a small command-line tool for converting between
  SAM/BAM/CRAM, primarily for benchmarking and exercising profiles.
- `CramStats`: dump per-block compression statistics from a CRAM file
  for debugging and ratio analysis.

Tests and CI
------------
- New hts-specs CRAM 3.0/3.1 compliance tests covering decode, index
  query, and round-trip behaviour.
- FQZComp round-trip tests using the hts-specs quality data files.
- CRAI index query correctness tests; codec roundtrip property tests.
- Split CRAM31 fidelity tests into per-profile classes for parallel
  execution.
- Reduced memory pressure in the unit tests to eliminate intermittent
  OOM failures on CI.
- Sped up several long-running CRAM tests by caching test data,
  reducing slice-size matrices, and downsampling fixtures.
- Misc: fixed a thread-safety bug in VariantContextTestProvider that
  was producing non-deterministic test counts.

CHANGELOG / README updated for the 5.0.0 release.

Co-authored-by: Chris Norman <cnorman@broadinstitute.org>

Author: tfenne

CHANGELOG.md

@@ -10,6 +10,48 @@ early infrastructure for a plugin-based codec framework and resource bundles.
 
 ---
 
+## 5.0.0
+
+Adds **CRAM 3.1 write support** to htsjdk.  This is the culmination of the read-side codec work
+in 4.2.0 and the reader wiring in 4.3.0: htsjdk can now produce CRAM 3.1 files that are
+interoperable with samtools/htslib.
+
+### CRAM 3.1 Write Support
+
+- Enable CRAM 3.1 writing with all spec codecs: rANS Nx16, adaptive arithmetic Range coder, FQZComp, Name Tokenisation, and STRIPE
+- Add configurable compression profiles (FAST, NORMAL, SMALL, ARCHIVE) with trial compression for automatic codec selection
+- Implement `TrialCompressor` to replace ad-hoc triple-compression for tags and align trial candidates with htslib
+- Add `GzipCodec` for direct Deflater/Inflater GZIP compression, wired into CRAM as a codec option
+- Strip NM/MD tags on CRAM encode and regenerate on decode, matching htslib behavior
+- Implement attached (same-slice) mate pair resolution
+- Align DataSeries content IDs with htslib for cross-implementation debugging
+- Remove content digest tags (BD/SD/B5/S5/B1/S1) from CRAM slice headers, matching htslib/samtools behavior. These are optional per the spec and were expensive to compute. Block-level CRC32 (required by CRAM 3.0+) provides data integrity. This is technically a breaking change but has zero practical impact since no known tools consume these tags.
+- Default CRAM version for writing is now 3.1 (was 3.0)
+- Add `CramConverter` command-line tool for testing and benchmarking CRAM write profiles
+
+### Codec and Compression Optimizations
+
+- Refactor and optimize all rANS codecs: byte-array API, backwards-write encoding, and general simplifications
+- Optimize Name Tokeniser encoder: replace regex with hand-written parser; add per-type flags, STRIPE support, stream deduplication, and all-MATCH elimination
+- Optimize FQZComp, Range coder, and rANS encoder hot paths
+- Tune NORMAL profile codec assignments based on empirical compression testing
+
+### Performance
+
+- Replace `ByteArrayInputStream`/`ByteArrayOutputStream` with unsynchronized `CRAMByteReader`/`CRAMByteWriter` to eliminate synchronization overhead
+- Fuse read base restoration, CIGAR building, and NM/MD computation into a single pass during decode
+- Cache tag key metadata to eliminate per-record `String` allocation during CRAM decode
+- Pool `RANSNx16Decode` instances in the Name Tokeniser
+- Optimize BAM nibble-to-ASCII base decoding with a bulk lookup table
+
+### Testing and Infrastructure
+
+- Split CRAM 3.1 fidelity tests into per-profile classes for parallel execution
+- Reduce memory pressure in unit tests to eliminate OOM failures
+- Fix thread-safety bug in `VariantContextTestProvider` causing non-deterministic test counts
+
+---
+
 ## 4.3.0 (2025-05-09)
 
 Completes CRAM 3.1 read support by wiring the codec implementations (added in 4.2.0) into

README.md

@@ -11,6 +11,8 @@ manipulating HTS data.
 
 > **NOTE: _HTSJDK has only partial support for the latest Variant Call Format Specification.  VCFv4.3 can be read but not written, VCFv4.4 can be read in lenient mode only, and there is no support for BCFv2.2._**
 
+> **NOTE: _HTSJDK now supports both reading and writing CRAM 3.1 files.  CRAM 3.1 write support includes all codecs defined in the specification (rANS Nx16, adaptive arithmetic Range coder, FQZComp, Name Tokenisation, and STRIPE), configurable compression profiles (FAST, NORMAL, SMALL, ARCHIVE), and trial compression for automatic codec selection.  Files produced by htsjdk are interoperable with samtools/htslib._**
+
 ### Documentation & Getting Help
 
 API documentation for all versions of HTSJDK since `1.128` are available through [javadoc.io](http://www.javadoc.io/doc/com.github.samtools/htsjdk).

build.gradle

@@ -99,7 +99,7 @@ tasks.withType(Test).configureEach { task ->
 
     // set heap size for the test JVM(s)
     task.minHeapSize = "1G"
-    task.maxHeapSize = "6G"
+    task.maxHeapSize = "12G"
 
     task.jvmArgs '-Djava.awt.headless=true'  //this prevents awt from displaying a java icon while the tests are running
 
@@ -160,7 +160,7 @@ test {
            excludeGroups "slow", "broken", "defaultReference", "optimistic_vcf_4_4", "ftp", "http", "sra", "ena", "htsget", "unix"
        }
        parallel = "classes"
-       threadCount = 2 * Runtime.runtime.availableProcessors()
+       threadCount = Runtime.runtime.availableProcessors()
    }
 } dependsOn testWithDefaultReference, testWithOptimisticVCF4_4
 
@@ -213,6 +213,12 @@ spotbugsTest {
     }
 }
 
+// Fat JAR with all dependencies for standalone CLI tools (CramConverter, CramComparison, etc.)
+shadowJar {
+    archiveClassifier.set('all')
+    mergeServiceFiles()
+}
+
 publishing {
     publications {
         htsjdk(MavenPublication) {

src/main/java/htsjdk/samtools/BamConverter.java

@@ -0,0 +1,122 @@
+package htsjdk.samtools;
+
+import java.io.File;
+
+/**
+ * Simple command-line tool for reading and optionally converting BAM files, primarily
+ * for experimenting with BAM read-path profiling.
+ *
+ * <p>Usage:
+ * <pre>
+ *   java -cp htsjdk.jar htsjdk.samtools.BamConverter input.bam [output.bam]
+ * </pre>
+ *
+ * <p>If no output is specified, records are read and iterated but not written.
+ */
+public class BamConverter {
+
+    private static final String USAGE = String.join("\n",
+            "Usage: BamConverter <input> [output]",
+            "",
+            "Read and optionally convert a BAM file.",
+            "",
+            "Arguments:",
+            "  input              Input BAM file",
+            "  output             Optional output BAM file (omit to read-only)"
+    );
+
+    /**
+     * Entry point. Parses command-line arguments and performs the read/conversion.
+     *
+     * @param args command-line arguments (see USAGE for details)
+     */
+    public static void main(final String[] args) {
+        if (hasFlag(args, "--help") || hasFlag(args, "-h")) {
+            System.out.println(USAGE);
+            System.exit(0);
+        }
+        if (args.length < 1) {
+            System.err.println(USAGE);
+            System.exit(1);
+        }
+
+        final boolean eager = hasFlag(args, "--eager");
+        // Collect positional args (non-flag arguments)
+        final String[] positional = java.util.Arrays.stream(args)
+                .filter(a -> !a.startsWith("--"))
+                .toArray(String[]::new);
+        if (positional.length < 1) {
+            System.err.println(USAGE);
+            System.exit(1);
+        }
+        final String inputPath = positional[0];
+        final String outputPath = positional.length > 1 ? positional[1] : null;
+
+        if (outputPath != null) {
+            System.err.printf("Converting %s -> %s%s%n", inputPath, outputPath, eager ? " (eager decode)" : "");
+        } else {
+            System.err.printf("Reading %s (no output%s)%n", inputPath, eager ? ", eager decode" : "");
+        }
+
+        final SamReaderFactory readerFactory = SamReaderFactory.makeDefault()
+                .validationStringency(ValidationStringency.SILENT);
+
+        long count = 0;
+        final long startTime = System.currentTimeMillis();
+
+        try (final SamReader reader = readerFactory.open(new File(inputPath))) {
+            final SAMFileHeader header = reader.getFileHeader();
+
+            if (outputPath != null) {
+                final SAMFileWriterFactory writerFactory = new SAMFileWriterFactory();
+                try (final SAMFileWriter writer = writerFactory.makeBAMWriter(header, true, new File(outputPath).toPath())) {
+                    for (final SAMRecord record : reader) {
+                        if (eager) record.eagerDecode();
+                        writer.addAlignment(record);
+                        count++;
+                        if (count % 1_000_000 == 0) {
+                            System.err.printf("  ... %,d records%n", count);
+                        }
+                    }
+                }
+            } else {
+                for (final SAMRecord record : reader) {
+                    if (eager) record.eagerDecode();
+                    count++;
+                    if (count % 1_000_000 == 0) {
+                        System.err.printf("  ... %,d records%n", count);
+                    }
+                }
+            }
+        } catch (final Exception e) {
+            die("Error: " + e.getMessage());
+        }
+
+        final long elapsed = System.currentTimeMillis() - startTime;
+        final long inputSize = new File(inputPath).length();
+
+        if (outputPath != null) {
+            final long outputSize = new File(outputPath).length();
+            System.err.printf("Done. %,d records in %.1fs. Input: %,d bytes, Output: %,d bytes (%.1f%%)%n",
+                    count, elapsed / 1000.0, inputSize, outputSize,
+                    inputSize > 0 ? (100.0 * outputSize / inputSize) : 0);
+        } else {
+            System.err.printf("Done. %,d records in %.1fs. Input: %,d bytes%n",
+                    count, elapsed / 1000.0, inputSize);
+        }
+    }
+
+    private static boolean hasFlag(final String[] args, final String flag) {
+        for (final String arg : args) {
+            if (flag.equals(arg)) return true;
+        }
+        return false;
+    }
+
+    private static void die(final String message) {
+        System.err.println("ERROR: " + message);
+        System.err.println();
+        System.err.println(USAGE);
+        System.exit(1);
+    }
+}

src/main/java/htsjdk/samtools/CRAMContainerStreamWriter.java

@@ -2,7 +2,7 @@
 
 import htsjdk.samtools.cram.build.ContainerFactory;
 import htsjdk.samtools.cram.build.CramIO;
-import htsjdk.samtools.cram.common.CramVersions;
+import htsjdk.samtools.cram.common.CRAMVersion;
 import htsjdk.samtools.cram.ref.CRAMReferenceSource;
 import htsjdk.samtools.cram.structure.*;
 import htsjdk.samtools.util.RuntimeIOException;
@@ -19,6 +19,7 @@ public class CRAMContainerStreamWriter {
     private final SAMFileHeader samFileHeader;
     private final ContainerFactory containerFactory;
     private final CRAMIndexer cramIndexer;
+    private final CRAMVersion cramVersion;
 
     private long streamOffset = 0;
 
@@ -70,7 +71,7 @@ public CRAMContainerStreamWriter(
      * Create a CRAMContainerStreamWriter for writing SAM records into a series of CRAM
      * containers on output stream, with an optional index.
      *
-     * @param encodingStrategy encoding strategy values
+     * @param encodingStrategy encoding strategy values (includes CRAM version)
      * @param referenceSource reference cramReferenceSource
      * @param samFileHeader {@link SAMFileHeader} to be used. Sort order is determined by the sortOrder property of this arg.
      * @param outputStream where to write the CRAM stream.
@@ -88,6 +89,7 @@ public CRAMContainerStreamWriter(
         this.outputStream = outputStream;
         this.cramIndexer = indexer;
         this.outputStreamIdentifier = outputIdentifier;
+        this.cramVersion = encodingStrategy.getCramVersion();
         this.containerFactory = new ContainerFactory(samFileHeader, encodingStrategy, referenceSource);
     }
 
@@ -103,11 +105,11 @@ public void writeAlignment(final SAMRecord alignment) {
     }
 
     /**
-     * Write a CRAM file header and the previously provided SAM header to the stream.
+     * Write a CRAM file header and the provided SAM header to the stream.
+     * Retained for backward compatibility with external projects (disq, GATK).
      */
-    // TODO: retained for backward compatibility for disq in order to run GATK tests (remove before merging this branch)
     public void writeHeader(final SAMFileHeader requestedSAMFileHeader) {
-        final CramHeader cramHeader = new CramHeader(CramVersions.DEFAULT_CRAM_VERSION, outputStreamIdentifier);
+        final CramHeader cramHeader = new CramHeader(cramVersion, outputStreamIdentifier);
         streamOffset = CramIO.writeCramHeader(cramHeader, outputStream);
         streamOffset += Container.writeSAMFileHeaderContainer(cramHeader.getCRAMVersion(), requestedSAMFileHeader, outputStream);
     }
@@ -131,7 +133,7 @@ public void finish(final boolean writeEOFContainer) {
                 writeContainer(container);
             }
             if (writeEOFContainer) {
-                CramIO.writeCramEOF(CramVersions.DEFAULT_CRAM_VERSION, outputStream);
+                CramIO.writeCramEOF(cramVersion, outputStream);
             }
             outputStream.flush();
             if (cramIndexer != null) {
@@ -144,7 +146,7 @@ public void finish(final boolean writeEOFContainer) {
     }
 
     protected void writeContainer(final Container container) {
-        streamOffset += container.write(CramVersions.DEFAULT_CRAM_VERSION, outputStream);
+        streamOffset += container.write(cramVersion, outputStream);
         if (cramIndexer != null) {
             // using silent validation here because the reads have been through validation already or
             // they have been generated somehow through the htsjdk

src/main/java/htsjdk/samtools/CRAMFileReader.java

@@ -374,7 +374,11 @@ public CloseableIterator<SAMRecord> queryUnmapped() {
         try {
             seekableStream.seek(0);
             iterator = new CRAMIterator(seekableStream, referenceSource, validationStringency);
-            seekableStream.seek(startOfLastLinearBin >>> 16);
+            // When startOfLastLinearBin is -1, there are no mapped reads and the entire file is
+            // unmapped. In that case, iterate from the beginning (already at position 0).
+            if (startOfLastLinearBin != -1) {
+                seekableStream.seek(startOfLastLinearBin >>> 16);
+            }
             boolean atAlignments;
             do {
                 atAlignments = iterator.advanceToAlignmentInContainer(SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX, SAMRecord.NO_ALIGNMENT_START);

src/main/java/htsjdk/samtools/CRAMIterator.java

@@ -115,7 +115,11 @@ private BAMIteratorFilter.FilteringIteratorState nextContainer() {
                     compressorCache,
                     getSAMFileHeader());
             samRecordIterator = samRecords.iterator();
-            return BAMIteratorFilter.FilteringIteratorState.MATCHES_FILTER;
+            // A container may match the query but produce no records (e.g. a container with
+            // only a compression header and no slices). Skip to the next container in that case.
+            return samRecords.isEmpty()
+                    ? BAMIteratorFilter.FilteringIteratorState.CONTINUE_ITERATION
+                    : BAMIteratorFilter.FilteringIteratorState.MATCHES_FILTER;
         } else {
             return BAMIteratorFilter.FilteringIteratorState.CONTINUE_ITERATION;
         }

src/main/java/htsjdk/samtools/SAMFileWriterFactory.java

@@ -60,6 +60,7 @@ public class SAMFileWriterFactory implements Cloneable {
     private SamFlagField samFlagFieldOutput = SamFlagField.NONE;
     private Integer maxRecordsInRam = null;
     private DeflaterFactory deflaterFactory = BlockCompressedOutputStream.getDefaultDeflaterFactory();
+    private CRAMEncodingStrategy cramEncodingStrategy = new CRAMEncodingStrategy();
 
     /** simple constructor */
     public SAMFileWriterFactory() {
@@ -76,6 +77,7 @@ public SAMFileWriterFactory( final SAMFileWriterFactory other) {
         this.tmpDir = other.tmpDir;
         this.compressionLevel = other.compressionLevel;
         this.maxRecordsInRam = other.maxRecordsInRam;
+        this.cramEncodingStrategy = other.cramEncodingStrategy;
     }
 
     @Override
@@ -241,6 +243,30 @@ public SAMFileWriterFactory setSamFlagFieldOutput(final SamFlagField samFlagFiel
         return this;
     }
 
+    /**
+     * Set the {@link CRAMEncodingStrategy} to use when creating CRAM writers. Controls the CRAM version,
+     * compression profile, and per-data-series codec selection.
+     *
+     * <p>The default strategy uses the {@link htsjdk.samtools.cram.structure.CRAMCompressionProfile#NORMAL} profile.
+     * To use a specific profile:
+     * <pre>
+     *   factory.setCRAMEncodingStrategy(CRAMCompressionProfile.ARCHIVE.toStrategy());
+     * </pre>
+     *
+     * @param cramEncodingStrategy the encoding strategy to use for CRAM output
+     * @return this factory for chaining
+     */
+    public SAMFileWriterFactory setCRAMEncodingStrategy(final CRAMEncodingStrategy cramEncodingStrategy) {
+        if (cramEncodingStrategy == null) throw new IllegalArgumentException("CRAM encoding strategy was null");
+        this.cramEncodingStrategy = cramEncodingStrategy;
+        return this;
+    }
+
+    /** @return the current CRAM encoding strategy */
+    public CRAMEncodingStrategy getCRAMEncodingStrategy() {
+        return cramEncodingStrategy;
+    }
+
     /**
      * Create a BAMFileWriter that is ready to receive SAMRecords.  Uses default compression level.
      *
@@ -524,8 +550,14 @@ public CRAMFileWriter makeCRAMWriter(final SAMFileHeader header, final OutputStr
      * @return CRAMFileWriter
      */
     public CRAMFileWriter makeCRAMWriter(final SAMFileHeader header, final OutputStream stream, final Path referenceFasta) {
-        // create the CRAMFileWriter directly without propagating factory settings
-        return new CRAMFileWriter(stream, new ReferenceSource(referenceFasta), header, null);
+        return new CRAMFileWriter(
+                cramEncodingStrategy,
+                stream,
+                null, // no index
+                true, // presorted
+                new ReferenceSource(referenceFasta),
+                header,
+                null);
     }
 
     /**
@@ -670,23 +702,17 @@ private CRAMFileWriter createCRAMWriterWithSettings(
 
         final Path md5Path = IOUtil.addExtension(outputFile, ".md5");
         final CRAMFileWriter writer = new CRAMFileWriter(
+                cramEncodingStrategy,
                 createMd5File ? new Md5CalculatingOutputStream(cramOS, md5Path) : cramOS,
                 indexOS,
                 presorted,
                 referenceSource,
                 header,
                 outputFile.toUri().toString());
-        setCRAMWriterDefaults(writer);
 
         return writer;
     }
 
-    // Set the default CRAM writer preservation parameters
-    private void setCRAMWriterDefaults(final CRAMFileWriter writer) {
-        //TODO: set encoding params
-        //writer.setEncodingParams(new CRAMEncodingStrategy());
-    }
-
     @Override
     public String toString() {
         return "SAMFileWriterFactory [createIndex=" + createIndex + ", createMd5File=" + createMd5File + ", useAsyncIo="