Yet Another Methylation Encoder

clilatest ↗

Bit-packing for DNA methylation. Arrays and whole genomes in one binary format family, where analysis is bitwise — and so stays fast across a thousandfold of scale, from 28K probes to 29M CpGs.

GitHub· KnowYourCG· SeSAMe· cinderplot· Cite ↗

Fetch — nothing ships with the binary

No reference data is bundled: no coordinates, no feature sets, no example methylomes. yame fetch opens a browser over the catalogue — species, then platform or genome build, then its files.

yame fetch

yame fetch   •   ~/.local/share/yame
        TARGET          KIND                                 TAG        IN STORE
  ▪ HUMAN                                                                       
❯ ▾ [*] hg38            genome                               mixed ▰▱▱▱▱    5/47
  ├──    ●  cpg_nocontig.cr                                  v2          28.9 MB
  ├── ▸ [ ] KYCG        sets                                 v2    ▱▱▱▱▱    0/32
  ├── ▸  ✓  data        sets                                 v3    ▰▰▰▰▰     5/5
  ────────────────────────────────────────
  hg38   genome build

Arrows move, / filters by name, source, collection or title, x takes a row (a folder takes everything under it) and f fetches. What the store already holds shows as present and cannot be re-taken.

Files land in a digest-verified store shared by every tool in the suite — $YAME_DATA_HOME, default ~/.local/share/yame — so what you fetch here is what KnowYourCG and SeSAMe find. Everything else on this page runs against it.

The store is laid out the way the tree reads. What you see under a name is what lands under that name on disk, so hg38/data in the browser is $YAME_DATA_HOME/hg38/data, and a genome's knowledgebase simply nests inside it. Which upstream repo published a directory is a column in fetch -l, not a path component — one name, whether you are browsing, fetching or looking at the disk.

What it carries are the row spaces everything else is measured in. Each length is distinct, so the count is the identity — which is what lets -R and -m take a name instead of a path.

row space        rows   row space        rows
─────────────────────────────────────────────
hg38       29,401,795   HM450         486,427
mm10       21,867,837   MM285         287,692
mm39       21,889,506   MSA           284,309
EPICv2        937,690   Mammal40       38,607
EPIC          866,553   HM27           27,722

Piped or redirected, fetch prints the catalogue as TSV instead, so a script never blocks on a keystroke. Naming a target, filtering it, dry runs and the rest: yame fetch -h.

Data conversion — three kinds of data, text in and out

Almost everything in practice is one of three kinds of data, and they land in two of the seven formats. One table covers how each is stored and what pack reads to build it.

the data                 format           on disk     pack  one line per CpG
────────────────────────────────────────────────────────────────────────────────
methylation, M/U counts  3  M/U counts    1-8 B/site  -fm   M ⇥ U
methylation, binarized   6  set+universe  2 bits/row  -fd   call ⇥ covered
a CpG set                6  set+universe  2 bits/row  -fd   in set ⇥ in universe

pack is the only way in. Input is one row per CpG already in the reference's order — nothing is sorted or matched for you. Columns are tab-separated, and under -fd the first column must be exactly 1 to count as in-set while the second is any nonzero int.

yame fetch -c hg38/data/human_hg38_test.cg
yame fetch hg38/cpg_nocontig.cr   # the coordinates, into the store
yame unpack -R hg38 human_hg38_test.cg | awk '$4!=2' | head -5

chr1  659296  659298  1
chr1  779099  779101  0
chr1  795737  795739  1
chr1  933194  933196  0
chr1  958160  958162  1

Format 6 prints one column by default; -f -1 gives set and universe separately, which is what pack -f6 reads back:

yame unpack -f -1 human_hg38_test.cg | awk '$1!="NA"' | head -3

1  1
0  1
1  1

yame unpack -f -1 human_hg38_test.cg > rt.txt
yame pack -f6 rt.txt rt.cg  # and back again

Verified byte-identical on the round trip.

Visualization — look at it without leaving the shell

Rows are samples, columns are CpGs or windows. The reference is inferred from the row count, so -r alone is enough.

yame fetch -c hg38/data/human_hg38_celltypes.cg
yame fetch hg38/cpg_nocontig.cr   # the coordinates, into the store
yame hprint -r chr20:30000000-30001000 human_hg38_celltypes.cg

                      chr20:30000011-30000862 (10 CpGs)
                      |.........
                      30000011  
Oligodendrocyte       █...██..░█
Pancreas-Beta         █.....░█..
Blood-NK              ....██░█░█
Blood-Monocytes       ....███.░█

-w caps the width — wider regions are window-averaged rather than truncated — and -R alone gives the whole-genome view.

Multi-sample indexing — many samples, one file

Records concatenate, so a cohort is one file. The .idx beside it maps sample names to offsets, which is what makes them addressable by name.

yame fetch -c hg38/data/human_hg38_celltypes.cg
yame info human_hg38_celltypes.cg | cut -f2-7

Sample           Ncol      Nrow  Format  UnitBytes  Keys
Oligodendrocyte     4  29401795       6          2    NA
Pancreas-Beta       4  29401795       6          2    NA
Blood-NK            4  29401795       6          2    NA
Blood-Monocytes     4  29401795       6          2    NA

yame chunk -s 500000 human_hg38_celltypes.cg chunks/  # fixed row blocks
yame subset -o two.cg human_hg38_celltypes.cg Blood-NK Blood-Monocytes
yame info human_hg38_celltypes.cg | cut -f2 | tail -n +2 > names.txt
yame split -s names.txt human_hg38_celltypes.cg out_  # one file each
yame index -s names.txt human_hg38_celltypes.cg  # rebuild the index

Sample               Nrow  Format
Blood-NK         29401795       6
Blood-Monocytes  29401795       6

info -1 reports the first record of each file — the cheap way to see what a directory holds:

File                              Sample                               Nrow  Format
human_hg38_40_celltypes_chr20.cg  GSM5652176_Adipocytes-Z000000T7    773477       3
human_hg38_celltypes.cg           Oligodendrocyte                  29401795       6
human_hg38_immune_mixture.cg      1                                29401795       3
human_hg38_test.cg                1                                29401795       6
human_hg38_test.truth.cg          1                                29401795       6

Compare & summarize — the analysis command

With no mask, summary describes the file. With one it gives the overlap table — universe, query, mask, overlap, log2 odds ratio, beta, depth. The mask is named, not spelled out: the query's row count says which row space to resolve it in.

yame fetch -c hg38/data/human_hg38_test.cg
yame fetch hg38/KYCG -g CGI   # the mask, into the store
yame summary -m CGI human_hg38_test.cg

QFile               Query  MFile            Mask     N_univ  N_query  N_mask  N_overlap  Log2OddsRatio   Beta  Depth
human_hg38_test.cg      1  CGI.20220904.cm  OpenSea   29000    23857   23534      21149           3.17  0.899     NA
human_hg38_test.cg      1  CGI.20220904.cm  Shelf     29000    23857    1224       1076           0.67  0.879     NA
human_hg38_test.cg      1  CGI.20220904.cm  Shore     29000    23857    2074       1263          -1.74  0.609     NA

A format 2 mask yields one row per state, so -m ChromHMM tests every chromatin state at once. -b picks the set in the browser instead.

yame fetch -c hg38/data/human_hg38_immune_mixture.cg
yame fetch -c hg38/data/human_hg38_40_celltypes_chr20.cg
yame binarize -t 0.5 -c 3 -o calls.cg human_hg38_immune_mixture.cg  # M/U -> set + universe
yame subset -o a.cg human_hg38_40_celltypes_chr20.cg GSM5652176_Adipocytes-Z000000T7
yame subset -o b.cg human_hg38_40_celltypes_chr20.cg GSM5652181_Saphenous-Vein-Endothel-Z000000RM
yame pairwise -H 1 -c 5 -d 0.2 a.cg b.cg  # differential, fmt3 -> fmt6
yame rowop -o binasum human_hg38_40_celltypes_chr20.cg pseudobulk.cg  # collapse

Sample    Nrow  Format
     1  773477       3

That last one is how pseudobulks are made: subset the cells of a cluster, pipe into rowop.

Perturb & mask — for benchmarking

Two different jobs. mask is deterministic and file-driven — it decides which sites are valid. dsample and perturb are stochastic and seed-driven — they decide how much signal is left.

yame fetch -c hg38/data/human_hg38_test.cg
yame dsample -N 5000 -s 42 -o ds.cg human_hg38_test.cg

Query  N_univ  N_query   Beta
    0    5000     4087  0.817

yame perturb -p 0.05 -s 7 -o noisy.cg human_hg38_test.cg

Query  N_univ  N_query   Beta
    1   29000    22953  0.791

Set -s. Without a fixed seed neither is reproducible. mask takes its mask as a positional argument, not -m, and both files must have the same row count.

yame fetch -c hg38/KYCG -g blacklist
yame mask -o clean.cg human_hg38_test.cg Blacklist.20220304.cm

Data model — row space, records, and the seven packings

A YAME file is values in order and nothing else. Row i means whatever row i of the reference means, so intersecting a methylome with a feature set is a bitwise AND rather than a genomic join — fast, and indifferent to how many rows there are.

cpg_nocontig.cr   the reference: every CpG in the genome, in order

  row               1       2       3       4       5       6
  chr1:           10,468  10,470  10,483  10,489  10,493  10,497

  .cg methylome     1       0       1       1       0       1
  .cm feature set   0       1       1       0       0       1
  ───────────────────────────────────────────────────────────────
  AND               0       0       1       0       0       1
                                    ▲
                     row 3 is chr1:10,483 in every file, so an
                     overlap is a bitwise AND, not a genomic join

Position is not stored; it is implied. A data file carries no coordinates at all — row i means whatever row i of the reference means. That is what makes the files small and the intersections free: overlapping a methylome with a feature set is a bitwise AND, not a genomic join. It is also why row counts have to match, and why a file's row count identifies its reference well enough for -R to be inferred.

Two files are comparable only if they share a row space, and since each has a distinct length, a row count is an identity. That is what lets -R and -m take a name instead of a path. An Infinium manifest and a whole-genome CpG set are the same kind of object here, differing only in length.

one .cg file = BGZF frames; any record reachable without the rest

  ┌────────┬─────┬────────┬───────────────────────────────┐
  │  sig   │ fmt │   n    │             data              │
  │  8 B   │ 1 B │  8 B   │  packed per the format code   │
  └────────┴─────┴────────┴───────────────────────────────┘
  └─────────────── one record = one sample ───────────────┘

  [ sample 1 ][ sample 2 ][ sample 3 ][ … ]
  concatenated; the .idx beside it stores where each begins

A record is a header and a blob. An 8-byte signature, the one-byte format code, the row count n, then the data packed according to that code. Records are concatenated, so a multi-sample file is just more of them, and the .idx beside it stores where each one begins — which is why yame subset can pull one sample out of hundreds without decompressing the others.

The three kinds above are the common cases. The rest exist because a row sometimes has to carry something else: a fraction, a label, a plain yes/no, or the coordinates themselves.

the data                 format           one line per CpG      on disk       NA        pack  what to know
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
methylation, M/U counts  3  M/U counts    M ⇥ U                 1-8 B/site    M=U=0     -fm   a tagged record; only 2^31+ counts lose bits
methylation, binarized   6  set+universe  call ⇥ covered        2 bits/row    00        -fd   the 2nd bit says the site was covered at all
a CpG set                6  set+universe  in set ⇥ in universe  2 bits/row    00        -fd   query and background travel in one file
a beta value             4  beta          a float, or NA        4 B/value     negative  -fn   NA runs collapse, so sparse is cheap
a plain yes/no           0  binary        0 or 1                1 bit/row     none      -fb   any first character but 0 reads as 1
a state label            2  states        the whole line        dict + runs   a label   -fs   keys stored once in a table, then runs
a small int              1  byte          one char per line     3 B/run       none      -fc   only the first character of a line
coordinates (.cr)        7  coordinates   BED: chrom ⇥ start    delta varint  n/a       -fr   only chrom and start are read

Format 5 is obsolete: still decodable, no longer packable.

Every format lays the same row space down at a different density, which is the whole reason there is a family rather than one format. Two of them are fixed width and cost the same per CpG whatever the data says; the rest widen with the values they carry, so a format 3 store scales with sequencing depth rather than with the row count. These layouts are the ones in src/format*.c.

one byte, and how many CpGs it holds

  fmt 0  binary        │1│0│0│1│1│0│1│0│   8 CpGs · 1 bit each
  fmt 6  set+universe  │11 │10 │00 │11 │   4 CpGs · 2 bits each
  fmt 1  small int     │      37       │   1 CpG  · one whole byte

  hg38 is 29,401,795 rows: 3.5 MB at one bit, 7.0 MB at two.

A fixed-width format is a packing density. These three spend the same number of bits on every CpG, so the file size follows from the row count alone: a presence/absence mask does not need a byte per CpG, and a small integer cannot live in a bit.

Format 5 is obsolete: still decodable, no longer packable.

Format 3 — M/U counts, tagged variable width

  inflated   c->unit bytes a site, default 8 (M and U 32 bits each);
             pack -u takes 1..8, and 1 leaves only 4 bits apiece
  on disk    the low two bits tag the record
             01   1 byte   [ M(3)  | U(3)  | 01 ]      M,U ≤ 6
             10   2 bytes  [ M(7)  | U(7)  | 10 ]      M,U ≤ 126
             11   8 bytes  [ M(31) | U(31) | 11 ]      M,U ≤ 2^31-1
             00   2 bytes  [ run length (14) | 00 ]    a run of M = U = 0
             beyond 2^31-1 fitMU right-shifts M and U together until both
             fit, so the ratio survives and the absolute depth does not

Format 6 — set + universe, two bits a row

  inflated   byte s[i>>2], offset = (i&3)*2           (n+3)>>2 bytes
             set bit      = (b >> offset)     & 1
             universe bit = (b >> (offset+1)) & 1
  codes      00 outside the universe · 10 in universe, 0 · 11 in universe, 1
  on disk    identical; "compression" is only a flag

Two bits is the interesting one: the universe bit separates not measured from measured and zero — a distinction one bit cannot make and every enrichment test needs, which is why a query and its background travel in one file.

two bits per CpG: the universe bit, then the set bit

  ┌──────────────┬──────────────┬──────────────┐
  │     0 0      │     1 0      │     1 1      │
  │      NA      │   value 0    │   value 1    │
  ├──────────────┼──────────────┼──────────────┤
  │   outside    │  measured,   │  measured,   │
  │ the universe │  not in set  │  in the set  │
  └──────────────┴──────────────┴──────────────┘

  The universe bit separates absent from present-and-zero — a
  distinction one bit cannot make, and every enrichment test needs.

Two bits, three meanings. A CpG that was never measured and a CpG measured as unmethylated are different facts; format 0 has to conflate them, format 6 does not. That is why a query and its background can travel in one file, and why summary can compute an odds ratio from a single record.

Format 0 — binary, one bit a row

  inflated   bit (i&7) of s[i>>3]                     (n+7)>>3 bytes
  on disk    identical; format 0 defines no coding of its own

Format 1 — one byte a row, run-length encoded

  inflated   s[i], any value 0-255
  on disk    3-byte records  [ value (1 B) | length (uint16 LE) ]
             length 1..32767; a longer run is split across records

Format 2 — categorical states, key table then runs

  layout     [ key section ][ '\0' ][ data section ]
  keys       "state0\0state1\0...\0", stored once; row value v means keys[v]
  inflated   data section is one integer a row, c->unit bytes
  on disk    [ unit ][ (value, length) ... ], unit in {1,2,3,8},
             length uint16 LE 1..65535

Format 4 — one float a row, NA runs collapsed

  inflated   float array, 4 bytes a row; NA is a negative float (-1.0)
  on disk    32-bit words, the MSB tags each
             1    [ 1 | run length (31) ]   that many consecutive NAs
             0    [ 0 | float bits (31) ]   the value, sign bit dropped

Format 7 — coordinates, delta varint by chromosome

  stream     [ chrom '\0' ][ delta ][ delta ] ... [ 0xff ][ chrom '\0' ] ...
  deltas     0bbbbbbb                   1 byte,  7-bit delta
             10xxxxxx yyyyyyyy          2 bytes, 14-bit delta
             11...... + 7 more          8 bytes, 62-bit delta, big-endian
  each delta is loc(i) - loc(i-1); 0xff ends a chromosome block