CRAM v4 idea tracker

[jkbonfield]

CRAM v4 isn't on the cards yet, but it's good to keep track of any ideas in one place so if it ever happens we don't forget something.

• Index format improvements: .crai index improvements #137
• Removal of CORE block?
• Cull some of the codecs (mainly those to do with core, eg beta, subexp, etc)
• Additional codecs, if not already adding by then in a 3.x sub-release. Eg FSE, Huff0 at the fast end, ZSTD in the middle and maybe fqzcomp for qualities and custom name-LZ at the slower end (see io_lib cram_modules branch).
• Read names present/absent via another route, eg zero length read name => generate. Currently it is done via a flag but this also forces detached mode and has other implications.
• Possibility of transforms prior to codec.
  • Eg nibble packing (put two qualities into one byte before applying rans0, rans1, huff0, etc).
  • String delta (prefix / suffix removal).
• Ability to store confidence values in original orientation. This improves order-1 compression ratio for qualities and is very quick to achieve.
• Support for multiple embedded references per slice?
• Permit rANS table to be in compression header and maybe implement an order-2 codec? It'll be much slower due to memory usage, but potentially with say 100 1k read slices per container we get size efficiency coupled to strong random access capability.
• Potential for conditional decoding with potential for duplicated entries. Eg h37 & h38.
  • Each read has a bit field to indicate under which reference the data should be emitted (both for most cases).
  • Possibility of bit field per data series too? Eg cigar may change, but quality string doesn't.
  • How to get working with read-pairing (read number delta)?

[mp15]

The first 4 items there could be done via a point release of 3 surely?

[jkbonfield]

I'm not so sure on index format changes or removal of features, but new codecs for sure.

[vadimzalunin]

Index format is independent of the file format, ideally it should have it's own versions.

New codecs are fine in v3.x as long as they are well justified.

Also there is the matter of priorities, for example embedded multirefs may have higher impact.

[jkbonfield]

Additional idea based on observation from RNAseq data.

Employ a special code to mean "same entry as last record". In particular we often see data (eg tophat) where, say, fwd read is aligned in two locations and rev read only matches one site. In order to permit valid PNEXT/RNEXT/TLEN fields the same record is duplicated twice in a row - once per fwd read.

A simple test of a block of quality values shows around 4% saving simply by embedding "duplicate last line" code, while keeping speed very high as it's trivial to implement. The same logic may apply to other fields, such as read-names, cigar strings, etc.

[jkbonfield]

More things to consider

• Support for verbatim CIGAR ops, permitting zero sized ops and distinguishing between =, X and M.
• NM / MD present bit field. This permits recreating these only for records that orgiinally had them and not for every record (or none).
• Some way of separating RNEXT, PNEXT and TLEN from each other. Right now they're all computed on-the-fly or all stored verbatim. We often find RNEXT and PNEXT correct and TLEN out by one, but then pay the cost of storing RNEXT and PNEXT verbatim. We also don't have a way of storing * 0 0 for these fields without verbatim storing them and marking the read as detached (which has implications elsewhere).

[jkbonfield]

Exploratory work based on a small chunk of Illumina NovaSeq data:

CRAM v3, 10k reads/slice: 23.5Mb, 5.9s encode, 1.6s decode
CRAM v4, 10k reads/slice: 19.7Mb, 7.9s encode, 2.0s decode

This is entirely down to two new codecs - a dedicated read name encoder (around 40% smaller than gzip and 30% smaller than lzma) and a rANS variant with bit-packing and RLE built-in. There were few tags in this data (it's from Isaac) so the names + qualities dominate. On older bwa data with 40 qualities I was looking at only 4% or so saving (vs 16% on this data). This is slightly slower, but note it is faster (at both encode and decode) and small than using either bzip2 or lzma so it offers us a better position on the pareto frontier for the size vs time tradeoff.

The rANS variant implementation poses an interesting question. How do we handle this in a general purpose manner?

Currently it's just a hack - the "order" byte (0 vs 1) now has two additional bits to indicate whether RLE was performed or whether bit-packing was present. Eg on a single block of NovaSeq quals:

268.4 MB/s enc, 505.4 MB/s dec	 15100000 bytes -> 713819 bytes (r4x16pr -o193)
182.6 MB/s enc, 318.9 MB/s dec	 15100000 bytes -> 768045 bytes (curr CRAMv3; rans4c -o1)
 24.5 MB/s enc,  33.8 MB/s dec   15100000 bytes -> 676298 bytes (bsc -tTp -e2 -m0)

(Bsc included for comparison with best reasonable-time general purpose tool, but it's still an order of magnitude more CPU. It perhaps indicates adaptive order-2 or higher coding is worth doing though.)

RLE is obvious, although some implementations are better than others. My version is based on mespotine method with a lookup table to indicate which symbols gain through RLE and which do not. I encode run-lengths and run-symbols as two distinct data sets, rans compressed and concatenated. In this regard it's perhaps similar to BYTE_ARRAY_LEN which has two encoders for lengths and bytes. So should RLE itself be a separate codec which can be combined with anything? That feels nice and generic, but may be slower (untested).

Bit-packing is basically applying a map, eg ACGT -> values 0,1,2,3 and packing 4 symbols into a byte with 2 bits each. It handles 0, 1, 2, 4 or 8 bits per symbol, with 0 being a special case to indicate "all the same value".

Again bit-packing could be its own codec with a single encoder for the output - like BYTE_ARRAY_STOP has one part meta data and one part encoder.

Then we could nest them: eg BIT_PACK(4, {"A",0,"C",1,"G",2,"T",3}, RLE(EXTERNAL(10), EXTERNAL(11)) meaning map ACGT to 0123, pack nibbles together, run-length the result and store run lengths in block 10 and run literals in block 11.

[jkbonfield]

Also of interest is a tweak to rans to indicate the N-way multiplexing of coders. Currently this is hard-coded at 4, but see the speed of implementations here which vary from 350Mb/s to 1400Mb/s decodes: https://github.com/jkbonfield/rans_static

A 32-way AVX2 accelerated rans can really shift, but it's inappropriate for most data sets as it has 8x as much overhead when flushing (128 bytes instead of 16, although there are smarter ways of flushing that may reduce it to maybe 70-80 bytes). We could just have multiple codecs with differing levels of unrolling, but more general is an N-way parameter with fixed allowable values (eg 1, 4 or 32 say) so we can pick the appropriate multiplexing for the size of the block being compressed.

We don't gain speed trying to massively unroll small data blocks but pay a price. Fundamentally the algorithm is the same in all cases and can be done in a single naive loop (see rNx16b.c), but permitting smarter SIMD methods if we have appropriate hardware.

[jkbonfield]

Adding to the ever growing idea list...

MAP function: converting listed integer values to codes 0, 1, 2, 3, ...

Eg for BF field we may currently define it as:

         BF =>         EXTERNAL {15}

But an analysis shows just 23 values used in an example file:

$ samtools view ~/scratch/data/9827_2#49.cram|head -10000|awk '{print $2}'|sort|uniq -c|sort -nr
   2418 163
   2415 83
   2361 99
   2358 147
    135 121
    133 73
     31 185
     19 137
     16 65
     14 81
     14 129
     14 113
     12 145
      9 97
      8 161
      8 1187
      8 1107
      7 177
      6 1171
      6 1123
      5 1161
      2 1145
      1 1097

Mapping those values to integers in descending frequency order permits better encoding. Currently HUFFMAN, BETA, GAMMA etc all write to CORE block, but could in theory bit-stream out to an external block instead. However we can just use EXTERNAL encoding once mapped and chuck it through a different entropy encoder. Eg in cram_dump output I'd expect to see something like:

         BF =>         MAP{23, {163,83,99,147,121,73,185,137,65,81,129,113,145,97,161,1187,1107,177,1171,1123,1161,1145,1097}, EXTERNAL {15}}

The effect of this is to pass all the ITF8 encoding inefficiencies from the external block to the compression header. Every time we emit 163 via EXTERNAL codec we are writing out two ITF8 bytes (0x80 0xa3) instead of one byte (0x00 after map), similarly for 147 (0x80 0x93 -> 0x3). This is an immediate size reduction, which can then also be further reduced by entropy encoders. A quick demo:

$ perl -e '%x=(163=>0,83=>1,99=>2,147=>3,121=>4,73=>5,185=>6,137=>7,65=>8,81=>9,129=>10,113=>11,145=>12,97=>13,161=>14,1187=>15,1107=>16,177=>17,1171=>18,1123=>19,1161=>20,1145=>21,1097=>22); while (<>) {chomp($_);print chr($x{$_})}' _BF_ |entropy16
len = 10000 bytes, entropy16  = 2867.470884, 2.293977 bits per byte
len = 10000 bytes, entropy8   = 2890.973451, 2.312779 bits per byte
len = 10000 bytes, entropy8o1 = 2861.055669, 2.288845 bits per byte

$ perl -e 'while (<>) {if ($_ >= 128) {print chr(0x80+($_>>8)),chr
($_&0xff)} else {print chr($_)}}' _BF_ | entropy16 
len = 14903 bytes, entropy16  = 3716.031005, 1.994783 bits per byte
len = 14903 bytes, entropy8   = 4592.193830, 2.465111 bits per byte
len = 14903 bytes, entropy8o1 = 2879.035860, 1.545480 bits per byte

In practice it's easy to achieve good entropy encoding for order-0 after mapping (we get 2946) and it's hard to overcome the order-1 frequency table cost when trying order-1 encoding of the itf8 values (we get 3224). So this tiny one-block demo is 8.6% saving of BF values, although they're a tiny portion of the overall file size. It's work on CF too, maybe some auxiliary tags, but probably little else. Therefore I don't have strong hopes for this one, but documenting it anyway just in case.

My thought is this also makes it easier to implement basic algorithms that work on byte streams at the external block level (eg RLE) that don't need to be more complex in order to handle streams of integers of unknown sizes. We can't tell if this is likely to be significant until it's tried out in practice.

[jkbonfield]

One more thought - ITF8 sucks for signed values! Eg reference ID is -2 (multi-ref via RI data series), -1 (unmapped) or >= 0 (reference N). -1 has to be stored of 0xff 0xff 0xff 0xff 0x0f. The final 0x0f ruins RLE encoding too (maybe that's why the original Java code set the unused bits to 1 instead of 0?).

A better solution, for the data series that have signed values, is to have an SITF8 encoding where the sign bit is moved from the top bit to the bottom bit, or zig-zag encoding (eg (n << 1) ^ (n >> 31) encode and (n >> 1) ^ -(n & 1) decode) to avoid having both +0 and -0 possibilities. This makes encoding of small negative numbers as efficient as encoding small positive numbers, instead of making all negatives take up 5 bytes. Note it's still beneficial to use the unsigned variant for something which can only be unsigned (eg number of records, length of span, etc).

[cmnbroad]

The v3 spec reserves 20 bytes in the header for a cramID, and suggests populating it with the filename as one possibility. The htsjdk implementation places the full pathname there (well, the first 20 bytes...). This causes the md5 of the file to vary depending on it's location when it was created. The htsjdk implementation can be changed, but we might want to also change the spec to be more prescriptive about what goes in here.

[vadimzalunin]

I think the purpose was sha1 digest of whatever, like the path. Vadim

…

On Thu, 11 Jan 2018, 11:07 Chris Norman, ***@***.***> wrote: The v3 spec reserves 20 bytes in the header for a cramID <https://github.com/samtools/hts-specs/blob/master/CRAMv3.tex#L402>, and suggests populating it with the filename as one possibility. The htsjdk implementation places the full pathname there (well, the first 20 bytes...). This causes the md5 of the file to vary depending on it's location when it was created. The htsjdk implementation can be changed, but we might want to also change the spec to be more prescriptive about what goes in here. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#144 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAoz5HUbcIEcp1NvoJlVtgdAxUnXfUBkks5tJjHHgaJpZM4IFXtD> .

[jkbonfield]

Sorry for the bump. Other random ideas that came up recently.

• Reference compression in sam_comp builds a model per locus and compares each read against this model instead of against the reference. This is dynamic and progresses as each successive read is added. Each locus has probability of A, C, G, T and *.

• Deez uses a similar idea but as a static model rather than as an adaptive model and with one call only (the "consensus"). Reads use the consensus to delta against and then the consensus is delta vs actual reference. Embedded ref in CRAM can already be used to do the former part, but not the latter part.

Merging both of these existing ideas together could be to have multiple embedded references per slice representing observed haplotypes and diffing reads aginst their closest haplotype. It's not a full statistical model, but does have the benefit of handling correlations. We'd need to use one of the standard and well tested clustering / compression approaches for describing how the multiple haplotypes compare against one another, with one canonical one (or their consensus) then comparing vs the reference.

This also has parallels with handling graph genomes, as that is simply a different approach for haplotype description. Reads are diffed against a graph node and graph nodes themselves are anchored on the reference.

[jkbonfield]

SlimGene uses reference based compression and, like CRAM, stores successive deltas for the position of each edit in the sequence. Unlike CRAM the positions are stored in order from 3' to 5' under the assumption that most short read data is high quality at the start and errors cluster towards the end. Ie the magnitude of the first_edit_pos is usually larger than the magnitude of (len - last_edit_pos).

I am unsure of how much difference this actually makes, but CRAM4 a protoype on NovaSeq data has the FP field consuming 6.3% of the total storage cost. Needs experimentation!