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. A folder is confirmed before it transfers — off a terminal that means -y, since a short name can reach a lot — but a name that picks out one file is its own confirmation. A filtered directory is still a directory, though, even when the filter leaves exactly one file, so the two -g lines below carry -y; with those, every fetch line on this page runs in a script as it stands. 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. With neither, every row is dumped as one glyph, which is how a whole reconstruction is compared against a whole truth: M/U counts and fractions render as H/M/L (-g for deciles), binary calls as 1/0/2.

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        NSample      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      mac70_mono30_2pow22              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 -y -g CGI hg38/KYCG   # 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
human_hg38_test.cg      1  CGI.20220904.cm  Island    29000    23857    2168        369          -5.10  0.170     NA

Read the last row first: CpG islands are depleted in this methylome (log2 odds −5.10) and what little overlaps is unmethylated (beta 0.170), against 0.899 in the open sea. That contrast is the whole point of the table, and it is one -m away.

Two spellings, two namespaces. fetch takes the BROWSER path, hg38/KYCG/CGI, because that is where the file goes in the store. -m takes the BARE name, CGI, and resolves it inside the query's own row space. The two are not interchangeable: -m hg38/KYCG/CGI fails, because the slashes make it look for a file of that literal name. The block above uses both, one line apart, and that is the asymmetry most likely to cost you an afternoon.

Many masks cost little more than one. The query is read once and every mask is measured against it 64 rows at a time, so a whole knowledgebase is one pass rather than one per set: the 1359-set TFBS collection against a 29.4M-row methylome takes 6 s where it used to take 193, with the same numbers and the same 45 MB. YAME_SUMMARY_PATH=1 reports on stderr how many masks took that path, since the output cannot tell you.

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 pairwise -S a.cg b.cg  # or score the agreement instead: one TSV line
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.

pairwise -S answers the other question about the same two samples — not where they differ but how much they agree — as one line of n mae rmse acc pearson over the sites both cover. Either side may be beta (format 4) rather than counts, since comparing a model’s output to its ground truth is inherently cross-format. Give it one sample and a multi-sample store and it scores every record against that sample, inflating the sample once while the store streams:

yame pairwise -S -1 GSM5652176_Adipocytes-Z000000T7 human_hg38_40_celltypes_chr20.cg predictions.cg

To score imputation rather than recall, blank the sites the model was shown before comparing — mask with the query as the mask leaves exactly the sites it had to guess:

yame mask truth.cg query.cg | yame pairwise -S -1 <cell> - predictions.cg

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
yame summary ds.cg          # dsample itself prints nothing

QFile  Query  MFile  Mask    N_univ  N_query  N_mask  N_overlap  Log2OddsRatio   Beta  Depth
ds.cg      0     NA  global    5000     4087    5000       4087             NA  0.817     NA

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

QFile     Query  MFile  Mask    N_univ  N_query  N_mask  N_overlap  Log2OddsRatio   Beta  Depth
noisy.cg      1     NA  global   29000    22953   29000      22953             NA  0.791     NA

Set -s. Without a fixed seed neither is reproducible, and the same seed reproduces a run only on the same build and platform — these counts come from Linux, and a macOS build of the same version draws different sites. mask takes its mask as a positional argument, not -m, and both files must have the same row count.

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

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

The four extensions. They are all the same container and differ only in what they hold, so any command reads any of them: .cg a methylome, one value per CpG; .cm a feature set or mask; .cr the coordinates of a row space; .cx the generic name, used where the content is not fixed — split writes it for that reason. Nothing enforces them; they are a convention for readers.

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

Full usage, the trust model for fetched data, and the format reference are in llms.txt; per-command options are in yame <command> -h. AGPL-3.0-or-later. Supported by NIH/NIGMS 5R35GM146978.