• BEDMatrix
  • Introduction to the BEDMatrix package
  • Limitations of the BEDMatrix package
• ALTREP
  • Introduction to the ALTREP framework
  • Example ALTREPs
• ALTREP Reimplementation of BEDMatrix
  • Step-by-Step
  • Demo
  • Conclusions

• A package to extract data from PLINK .bed files

• .bed files contain genotype calls at biallelic variants

  library(BEDMatrix)
  path <- system.file("extdata", "example.bed", package = "BEDMatrix")
  X <- BEDMatrix(path)
  dim(X) # matrix-like object

  ## [1]   50 1000

  X[1:5, 1:5]

  ##           snp0_A snp1_C snp2_G snp3_G snp4_G
  ## per0_per0      0      1      1      1      0
  ## per1_per1      1      1      1      1     NA
  ## per2_per2      1      0      0      2      0
  ## per3_per3      2      0      0      0      1
  ## per4_per4      0      1      0      0      0

• 0x6c 0x1b (magic number) 0x01 (mode: variant-major)
• Sequence of V blocks of N/4 (rounded up) bytes each, where V is the number of variants and N is the number of samples
• The low-order two bits of a block's first byte store the first sample's genotype code. The next two bits store the second sample's genotype code, and so on for the 3rd and 4th samples. The second byte stores genotype codes for the 5th-8th samples, the third byte stores codes for the 9th-12th, etc.
• The two-bit genotype codes have the following meanings:
  • 00: Homozygous for first allele in .bim file
  • 01: Missing genotype
  • 10: Heterozygous
  • 11: Homozygous for second allele in .bim file
• If N is not divisible by four, the extra high-order bits in the last byte of each block are always zero.

Supplementary files:

• .fam file contains sample information (including number of samples)
• .bim file contains variant information (including number of variants)

See https://www.cog-genomics.org/plink/1.9/formats#bed

• Currently implemented as an S4 class that contains an external pointer to a C++ class that contains a pointer to the memory-mapped .bed file and provides methods for extracting genotypes as integer vectors

• Access is read-only

  as.vector(class(X))

  ## [1] "BEDMatrix"

  slotNames(X)

  ## [1] "xptr"   "dims"   "dnames" "path"

  methods(class = "BEDMatrix") # lots of methods need to be reimplemented

  ##  [1] as.matrix  [          dim        dimnames<- dimnames   initialize
  ##  [7] is.matrix  length     show       str       
  ## see '?methods' for accessing help and source code

Some common issues:

• Internal and primitive functions without methods dispatch will likely not work (hence BGData)

  tcrossprod(X)

  ## Error in tcrossprod(x, y): requires numeric/complex matrix/vector arguments

• Saving and reloading an object does not work

  save(X, file = "X.RData")
  load("X.RData")
  X[1]

  ## Error in extract_vector(x, i, ...): external pointer is not valid

• …
• Arithmetic sequences created by 1:n, seq_along, and the like now use compact internal representations via the ALTREP framework.
• …
• R has new serialization format (version 3) which supports custom serialization of ALTREP framework objects.
• …
• R 3.5.0 includes a framework that allows packages to provide alternate representations of basic R objects (ALTREP). The framework is still experimental and may undergo changes in future R releases as more experience is gained. For now, documentation is provided in https://svn.r-project.org/R/branches/ALTREP/ALTREP.html.

Source: https://cran.r-project.org/doc/manuals/r-release/NEWS.html

• R objects: vectors of various modes, functions, environments, language objects, etc.
• For the sake of this talk, we have to look at implementation details:
  • In C code, R objects are represented by the types SEXP and VECSEXP, which are pointers to a structure with typedef SEXPREC or VECTOR_SEXPREC

Example: x <- 1:10

• In R code: integer vector
• In C code: VECTOR_SEXPREC of type INTSXP

Basic representation of x <- 1:10 in memory on a 64-bit machine:

• As a sequence: 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
  • With data in memory (like the basic representation)
  • With data from somewhere else (file, network, database, …)
• As three numbers: from = 1, to = 10, by = 1 with evaluation on demand
• …
• This becomes more important as data gets bigger

writeBin(1:10, con = "sequence.bin") # write example data (imagine something bigger)

mmap <- function(filename, type = c("integer", "double"),
                 ptrOK = TRUE, wrtOK = FALSE, serOK = TRUE) {
    type <- match.arg(type)
    .Internal(mmap_file(filename, type, ptrOK, wrtOK, serOK))
}


x <- mmap("sequence.bin", type = "integer")
x

##  [1]  1  2  3  4  5  6  7  8  9 10

typeof(x) # regular integer vector

## [1] "integer"

class(x) # no method magic

## [1] "integer"

x <- 1:1e12 # huge sequence
format(object.size(x), units = "GB") # seemingly too large to allocate

## [1] "7450.6 Gb"

x[7] # but can read from it

## [1] 7

typeof(x) # regular numeric vector

## [1] "double"

class(x) # no method magic

## [1] "numeric"

• New ALTREPs can be developed and distributed in the form of packages
• Currently, there is not much documentation about this mechanism
  • Documentation: https://svn.r-project.org/R/branches/ALTREP/ALTREP.html
  • Presentation: https://homepage.divms.uiowa.edu/~luke/talks/uiowa-2018.pdf
  • R source code: in particular src/include/R_ext/Altrep.h and src/main/altrep.c (R > 3.5.0)
  • ALTREP examples repo on GitHub: in particular Rpkg-simplemmap
• Could ALTREP solve some of the limitations and complexities of the current BEDMatrix approach?

• The mmap_file internal function is currently limited to files that store data as a sequence of 32-bit integers or doubles
• The .bed format is not structured that way
• Does that mean it's over? Let's try to create our own ALTREP called ALTBED and see how far we can go!
  • The code presented here is simplified and incomplete
  • Full code available on GitHub: https://github.com/agrueneberg/ALTBED

1. Create new R package (for example with devtools::create_package("ALTBED")) with some infrastructure to make use the R's C API
   • Create src/ directory
   • Create source file in src/ that contains an implementation of R_init_ALTBED()
   • Add useDynLib("ALTBED", .registration = TRUE) to NAMESPACE
2. Create ALTREP class in R_init_ALTBED()
3. Override ALTREP methods in R_init_ALTBED()
4. Instantiate ALTREP objects in registered routines
5. Add some glue code in R that calls these registered routines using .Call()

• Only ALTINTEGER, ALTREAL, and ALTSTRING can be subclassed by package developers at the moment
  • ALTINTEGER is an alternate representation of an integer vector (INTSXP)
  • ALTREAL is an alternate representation of a double vector (REALSXP)
  • ALTSTRING is an alternate representation of a character vector (STRSXP)
• In our case, ALTINTEGER is the best fit for representing genotype calls (0/1/2/missing)

#include <R.h>
#include <Rinternals.h>
#include <R_ext/Altrep.h>

static R_altrep_class_t altrep_class_bedmatrix;

void R_init_ALTBED(DllInfo *dll) {
 // ...
    R_altrep_class_t class = R_make_altinteger_class("BEDMatrix", "ALTBED", dll);
    altrep_class_bedmatrix = class;
 // ...
}

• Objects of this class are still SEXPs – with the alt bit set in sxpinfo
• At runtime, method dispatch will happen based on the presence of the alt bit (in a way, similar to S4 but at a much deeper level)

• ALTREP
  • Length: get number of data elements
  • Duplicate and DuplicateEX: duplicate SEXP
  • Coerce, Serialized_state, Unserialize, UnserializeEX, Inspect
• ALTVEC
  • Dataptr: get pointer to data elements
  • Dataptr_or_null: get pointer to data elements or NULL otherwise
  • Extract_subset: get data elements that match index
• ALTINTEGER
  • Elt: get ith element
  • Get_region: get n elements starting from i
  • Is_sorted, Min, Max, Sum, No_NA

Dataptr, Dataptr_or_null, Extract_subset, and Elt all extract data in different ways in different contexts:

• Dataptr is used with macros such as INTEGER() and returns a pointer to 32-bit ints in memory
  • In our case, we cannot meaningfully return a pointer to the memory-mapped file because of the structure of a .bed file (same reason why we cannot use mmap_file())
    • Technically, it appears that we could allocate an integer array and return to a pointer to that, but that defeats the purpose
• Dataptr_or_null returns a data pointer if possible, or NULL otherwise (used to test if a data pointer can be returned)
  • In our case, this should always return NULL to fall back to a different mechanism
• Extract_subset is used with subsetting operators such as [
  • In our case, this is possible to implement (oh joy: it even converts the indices)
• Elt is used for extracting a single data element at position k
  • In our case, this is also possible to implement (a simpler case of Extract_subset)
• Get_region is used for extracting a sequence of elements (in case it is more efficient to extract one by one)
  • Default uses Elt so we do not need to implement this

• Duplicate necessary to maintain the "call by value" illusion
  • By default, this serializes the data pointer
  • In our case, we cannot expose the data pointer, but the object cannot be modified anyway, so we should be OK with returning the original object (not sure about this one)
• Serialized_state is called on save() (currently, version = 3 needs to be specified with save() for ALTREP support)
  • In our case, it should return everything needed to reconstruct the object
• Unserialize is called on load()
  • In our case, it should reconstruct the object based on the serialized state
• Returning NULL from a method typically indicates that there is no special behavior and that the same behavior as for the basic representation should be used

static R_altrep_class_t altrep_class_bedmatrix;

void R_init_ALTBED(DllInfo *dll) {
 // ...
 // create class
    R_altrep_class_t class = R_make_altinteger_class("BEDMatrix", "ALTBED", dll);
    altrep_class_bedmatrix = class;
 // override ALTREP methods
    R_set_altrep_Length_method(class, bedmatrix_Length);
    R_set_altrep_Duplicate_method(class, bedmatrix_Duplicate);
 // override ALTVEC methods
    R_set_altvec_Dataptr_method(class, bedmatrix_Dataptr);
    R_set_altvec_Dataptr_or_null_method(class, bedmatrix_Dataptr_or_null);
    R_set_altvec_Extract_subset_method(class, bedmatrix_Extract_subset);
 // override ALTINTEGER methods
    R_set_altinteger_Elt_method(class, bedmatrix_Elt);
 // ...
}

static SEXP make_bedmatrix(SEXP path, SEXP nrows, SEXP ncols) {
 // ...
 // mmap .bed file
 // ...
 // Create ALTREP object
    SEXP state = PROTECT(make_bedmatrix_state(...));
    SEXP eptr = PROTECT(make_bedmatrix_eptr(...));
    SEXP altrep = PROTECT(R_new_altrep(altrep_class_bedmatrix, eptr, state));
 // Turn it into a matrix
    SEXP dim = PROTECT(Rf_allocVector(INTSXP, 2));
    INTEGER(dim)[0] = Rf_asInteger(nrows);
    INTEGER(dim)[1] = Rf_asInteger(ncols);
    Rf_setAttrib(altrep, R_DimSymbol, dim);
    UNPROTECT(4); // state, eptr, altrep, dim
    return altrep;
}

SEXP C_map(SEXP path, SEXP nrows, SEXP ncols) {
    return make_bedmatrix(path, nrows, ncols);
}

static const R_CallMethodDef CallEntries[] = {
    {"C_map", (DL_FUNC) &C_map, 3},
    {NULL, NULL, 0}
};

void R_init_ALTBED(DllInfo *dll) {
 // ...
    R_registerRoutines(dll, NULL, CallEntries, NULL, NULL);
 // ...
}

map <- function(path) {
    # Determine n from .fam file and p from .bim file ...
    .Call(C_map, path, n, p)
}

X <- map("inst/extdata/example.bed")

Pros:

• Proper serialization and unserialization support
• Some functions such as tcrossprod work, but load data into memory (chunked processing is still needed)
  • This may improve greatly as ALTREP matures (e.g., do_summary uses Get_region)
• Easier to implement subsetting (no need for crochet)

Cons:

• Objects cannot be printed (at least in our case because we don't allow Dataptr)
• Probably more bugs because it's all written in C :-)