# Using bigmemory in R for Genomic Matrices Too Large for RAM Genotype matrices outgrow RAM quietly. You run a QC script happily on a chromosome, then the same code on a whole-genome file dies with an allocation error. Population biobanks now routinely hold matrices of millions of variants measured on hundreds of thousands of individuals [1], and even a modest 50,000 marker array typed on 500,000 animals fills 200 gigabytes at double precision. Script-based environments such as R hit these walls hard, a limitation serious enough that the builders of BESSiE, a tool for genomic best linear unbiased prediction (GBLUP) and Bayesian [Markov chain Monte Carlo](/blog/data-analysis/markov-chain-monte-carlo-definition-examples) (MCMC) analyses of large-scale marker data, stepped outside plain R data structures entirely [2]. The bigmemory package lets you stay in R. It stores a matrix in a file on disk, maps that file into your R process, and hands you an object you index like an ordinary matrix while the operating system streams only the pieces you touch. This guide walks the full workflow on a genotype dosage matrix, where each cell holds the count of alternate alleles carried, coded 0, 1, or 2 for hard calls. You will create a file-backed matrix, reopen it in a fresh session, compute allele frequencies and quality metrics without loading the data, filter variants, score samples in bounded memory, run parallel workers against one shared matrix, and fit penalized regression models straight from disk. You need R 4.x or later, about 2 GB of free disk for the demo, and three CRAN packages: bigmemory, biganalytics, and biglasso. Everything runs on a laptop. The worked example simulates a small dosage matrix so you can verify every command, and each step notes what changes at real cohort scale. Large-matrix handling is one of the real dividing lines in R for bioinformatics, and the pattern here transfers to expression matrices, methylation arrays, and mutation count matrices as well. ## Why Genomic Matrices Outgrow RAM in R R keeps every object in main memory, and several common operations copy objects rather than view them. Reading a dosage CSV with read.csv builds a data.frame with character columns for identifiers and double columns for values. Coercing to a matrix copies the numeric data again. Building a model matrix copies it a third time. A matrix that occupies 100 GB on disk can demand 300 to 400 GB of transient RAM before a single model is fitted. The arithmetic of genomic matrices makes this worse because both dimensions are large: - A 50,000 SNP array on 500,000 animals holds 25 billion cells. At 1 byte per stored genotype that is 25 GB, at 8 bytes per double it is 200 GB. - Whole-genome sequence dosages at one million variants on the same cohort reach 500 billion cells, 500 GB even at 1 byte per cell. - A genomic relationship matrix, the square matrix of pairwise relatedness used in GBLUP, grows with the square of sample size. For 250,000 samples it needs 250,000 × 250,000 × 8 bytes, about 500 GB, before you compute a single element. [Genomic prediction](/knowledge/bioinformatics/genomic-prediction-in-livestock-a-decision-framework-for-breeders) pipelines feel this early. GBLUP on 50k or high-density SNP arrays requires that square relationship matrix [3], and Bayesian variable-selection models such as BayesA, BayesB, BayesC, and BayesR iterate over the marker matrix many times [2]. Representation choices matter more than hardware in this regime. The snpnet work on biobank-scale regression encodes each genotype in two bits rather than eight bytes, a 32-fold memory reduction that makes regularized regression on millions of variants feasible on ordinary hardware [1]. bigmemory attacks the problem at the storage layer. Instead of deciding what fits in RAM, it moves the matrix to disk and lets the operating system manage which parts stay resident. ## How bigmemory Works Under the Hood A big.matrix is a small R object holding a pointer to a C++ data array. Two modes exist. An in-memory shared big.matrix lives in RAM but can be shared between R processes on one machine. A file-backed big.matrix lives in a backing file on disk, and that is the mode you want for genomic matrices larger than RAM. Three files tell the story: - The backing file, for example geno.bin, is a flat binary array of nrow × ncol cells at the chosen element size. A 2,000 by 5,000 char matrix produces exactly 10,000,000 bytes. - The descriptor file, for example geno.desc, is a tiny metadata record holding dimensions, element type, and the backing file location. You share the descriptor, never the data. - Your R session calls attach.big.matrix on the descriptor and receives a big.matrix object pointing at the same backing. The operating system treats the backing as a memory-mapped file. Pages of the file enter RAM on demand, and the page cache, the kernel's pool of recently used disk pages, evicts them under pressure. R never needs the whole matrix resident. Reading a column touches consecutive pages, which disks handle efficiently. Storage is column-major, matching base R. Columns are contiguous on disk, so column statistics and column-blocked operations stream smoothly, while row-wise access on a wide matrix hops across the file. This single fact drives most of the workflow design advice below. Element types trade precision for space. The char type stores one byte per cell, enough for dosages coded 0, 1, 2. The double type stores eight bytes and is required by some modeling packages. The table in Step 1 lists all four. bigmemory is deliberately low level. It provides storage, indexing, and basic statistics, with no compression, no sparse structures, and no genotype-aware operations. Dense file backing is the wrong substrate for extremely sparse data such as single-cell count matrices, where sparse data structures cut memory overhead dramatically, the redesign that let CoGAPS 3 scale Bayesian matrix factorization to roughly 1,000-fold more cells [4]. For dense genotype panels, dense file backing is usually right. ## Setting Up Your R Environment ```r install.packages(c("bigmemory", "biganalytics", "biglasso")) ``` bigmemory provides the core classes, biganalytics adds column statistics and the bigglm regression fitter, and biglasso adds penalized regression against file-backed design matrices. Companion packages you may meet in older workflows include bigtabulate for cross-tabulation, bigalgebra for matrix algebra, and bigmemory.sri for shared resource management across sessions. bigmemory is distributed through CRAN rather than Bioconductor. If your pipeline is Bioconductor-centric, most Bioconductor tutorial material for on-disk matrices centers on DelayedArray with an HDF5 backend through the HDF5Array package, installed with BiocManager. Both approaches solve the same problem, and the choice usually follows the rest of your stack. The steps below use CRAN packages only. ## Step-by-Step: bigmemory in R for a Large Genomic Matrix The demo builds a 2,000 sample by 5,000 SNP dosage matrix. Every command scales to real cohorts by changing numbers and block sizes, not structure. ### Step 1. Size the Matrix and Choose a Storage Type Compute the backing size before you commit: ``` backing bytes = nrow × ncol × bytes per cell ``` For the demo, 2,000 × 5,000 × 1 byte = 10,000,000 bytes, about 10 MB. For 500,000 samples and 50,000 SNPs as char, 25 GB. As double, 200 GB. Pick the smallest type that represents your data exactly. | Type | Bytes per cell | Represents | Typical genomic use | Watch out for | |--|--|--|--|--| | char | 1 | Integers from -128 to 127 | Dosage calls 0, 1, 2 with 3 reserved for missing | Decimals truncate silently | | short | 2 | Integers from about -32,768 to 32,767 | Scaled dosages, for example dosage × 10,000 | Same truncation caution | | integer | 4 | Integers up to about 2.1 billion | Counts, scaled values, indexes | Heavier than needed for calls | | double | 8 | Full double precision | Continuous dosages, expression values, model frames | Four to eight times the disk footprint | Two block-sizing formulas cover most workflows. For column blocks over n rows, columns per block ≈ RAM budget in bytes ÷ (n × bytes per cell). With 250,000 rows of char and a 2 GB budget, about 2.1 billion bytes, that is roughly 8,000 columns per block. For row blocks over p columns, swap the roles. ### Step 2. Create the File-Backed Matrix Simulate a dosage matrix and write it to a file-backed big.matrix in one call: ```r library(bigmemory) set.seed(42) n <- 2000 # samples p <- 5000 # SNPs geno <- matrix( sample(0:2, n * p, replace = TRUE, prob = c(0.45, 0.35, 0.20)), nrow = n, ncol = p ) dir.create("bigdata", showWarnings = FALSE) bm <- as.big.matrix( geno, type = "char", backingfile = "geno_demo.bin", descriptorfile = "geno_demo.desc", dirpath = "bigdata/" ) dim(bm) bm[1:4, 1:6] file.size("bigdata/geno_demo.bin") ``` Representative output, your genotype values will differ with the random draw: ``` [2] 2000 5000 [,1] [,2] [,3] [,4] [,5] [,6] [1,] 0 2 1 2 0 1 [2,] 1 1 0 1 2 2 [3,] 2 0 1 0 1 0 [4,] 1 2 2 1 0 1 [2] 10000000 ``` The parameters that matter: type sets the element size, backingfile names the binary array, descriptorfile names the metadata record, and dirpath places both in one folder. as.big.matrix copies from an in-memory source, which is fine for the demo but impossible for matrices larger than RAM. For real data, stream straight from a delimited file with read.big.matrix so the full matrix never exists in RAM: ```r bm <- read.big.matrix( filename = "geno_dosages.csv", sep = ",", header = FALSE, type = "char", backingfile = "geno_dosages.bin", descriptorfile = "geno_dosages.desc", dirpath = "bigdata/" ) ``` read.big.matrix parses the text file in chunks and writes the backing as it goes, so peak RAM stays small regardless of file size. Text parsing is slow at terabyte scale, so prefer binary sources when you have them. If your export carries an identifier column, drop it before reading, and keep sample identifiers in a separate key file with one ID per row in file order. Character identifiers inside the matrix force a wider type, and the key file habit survives every later export. ### Step 3. Reattach the Matrix in a Fresh Session Close R, reopen it, and attach through the descriptor: ```r library(bigmemory) bm <- attach.big.matrix("bigdata/geno_demo.desc") dim(bm) bm[1:3, 1:3] ``` ``` [2] 2000 5000 [,1] [,2] [,3] [1,] 0 2 1 [2,] 1 1 0 [3,] 2 0 1 ``` The descriptor also travels as an R object: ```r saveRDS(describe(bm), "bigdata/geno_demo_descriptor.rds") # any later session bm <- attach.big.matrix(readRDS("bigdata/geno_demo_descriptor.rds")) ``` Always distribute the descriptor with the backing file as a pair. Saving the big.matrix object itself with saveRDS preserves a pointer that means nothing in a new session. ### Step 4. Run Column-Level Quality Control Allele frequency per SNP falls out of a column mean. For dosages coded 0, 1, 2, the alternate allele frequency is the column mean divided by 2: ```r library(biganalytics) afs <- big.colMeans(bm) / 2 summary(afs) ``` ``` Min. 1st Qu. Median Mean 3rd Qu. Max. 0.3421 0.3691 0.3749 0.3750 0.3809 0.4087 ``` The simulated allele frequency sits near 0.375 by construction, which is exactly the kind of closed-form check you should run against array expectations on real data. Check genotype coding on a single column. Extracting one column of 2,000 cells costs almost nothing: ```r table(bm[, 1]) ``` ``` 0 1 2 902 687 411 ``` Counts should sum to the number of samples, and unexpected codes such as negative values or stray 9s show up here immediately. To locate rows carrying a missing code, mwhich scans a column without extracting it: ```r missing_rows <- mwhich(bm, cols = 1, values = 3, ops = "eq") head(missing_rows) ``` The demo data carries no 3 codes, so the result is an empty integer vector. On real data this returns row indices worth inspecting. Scripted, standardized QC scales better than eyeballing. The GWAS community formalized this early, with GWAtoolbox checking formatting, quality statistics, and cross-study consistency across large numbers of meta-analysis files in minutes [5]. Apply the same discipline to your matrix: assert the coding, the missing rate, and the allele frequency range in code. ### Step 5. Filter Variants Without Loading the Matrix Filtering by minor allele frequency (MAF), the smaller of the two allele frequencies at a variant, produces a smaller matrix. The trap is writing bm[, keep], which extracts the surviving columns into RAM as a base R matrix. Instead, allocate a second file-backed matrix and copy in column blocks: ```r maf <- pmin(afs, 1 - afs) keep <- which(maf > 0.05) p_keep <- length(keep) bm_f <- big.matrix( nrow = nrow(bm), ncol = p_keep, type = "char", backingfile = "geno_maf.bin", descriptorfile = "geno_maf.desc", dirpath = "bigdata/" ) cols_per_block <- 500 for (start in seq(1, p_keep, by = cols_per_block)) { end <- min(start + cols_per_block - 1, p_keep) bm_f[, start:end] <- bm[, keep[start:end]] } dim(bm_f) ``` ``` [2] 2000 5000 ``` Every SNP passes the 0.05 threshold in the simulated data, so the filtered matrix matches the original. The mechanics are the point. On real sequencing data this step commonly drops a large fraction of rare variants, and the block size is your RAM dial. At 250,000 rows and char storage, 8,000 columns per block keeps each extraction near 2 GB, from the formula in Step 1. ### Step 6. Score Every Sample With Chunked Matrix Products A polygenic score is a weighted sum of a sample's dosages across SNPs. The weight vector is tiny, so the only constraint is streaming the matrix in column blocks and accumulating: ```r set.seed(7) beta <- rnorm(ncol(bm_f)) prs <- numeric(nrow(bm_f)) cols_per_block <- 1000 for (start in seq(1, ncol(bm_f), by = cols_per_block)) { end <- min(start + cols_per_block - 1, ncol(bm_f)) block <- bm_f[, start:end] # 2000 x 1000, about 8 MB prs <- prs + as.numeric(block %*% beta[start:end]) } head(prs) ``` ``` [2] 33.41 -81.27 12.09 -47.85 64.52 -20.93 ``` Values are representative. Production scores center genotypes per SNP by subtracting twice the allele frequency before applying weights, and take weights from a published or fitted source rather than random numbers. The pattern generalizes to any column-wise statistic: extract a block, reduce it, accumulate the result, discard the block. Peak RAM stays at the block size no matter how many samples sit on disk. ### Step 7. Analyze in Parallel Against One Shared Matrix Because the matrix lives in one backing file, worker processes attach to the same descriptor and read disjoint column ranges with no duplication: ```r library(foreach) library(doParallel) desc <- describe(bm_f) starts <- seq(1, ncol(bm_f), by = 1000) chunks <- lapply(starts, function(s) s:min(s + 999, ncol(bm_f))) cl <- makeCluster(4) registerDoParallel(cl) means_par <- foreach( chunk = chunks, .combine = "c", .packages = "bigmemory" ) %dopar% { bm_w <- attach.big.matrix(desc) big.colMeans(bm_w, cols = chunk) } stopCluster(cl) all.equal(means_par, big.colMeans(bm_f)) ``` ``` [2] TRUE ``` Each worker pays one cheap attach, reads its slice at disk speed, and returns a small vector. The all.equal check confirms the parallel result matches a serial pass, a habit worth keeping for every parallel rewrite. Use makeCluster rather than forked parallelism when the parent process already holds large objects. Forked sessions on macOS can end abruptly under memory pressure rather than throwing a catchable error, while PSOCK workers, separate R sessions connected over sockets, start clean. ### Step 8. Fit Models Directly From Disk Two fitters accept file-backed matrices. bigglm from biganalytics pushes rows through a generalized linear model in chunks, so phenotype and covariate columns travel with the genotypes. Build a small demo matrix with named columns: ```r y <- rnorm(n) age <- sample(20:80, n, replace = TRUE) sex <- sample(0:1, n, replace = TRUE) model_df <- data.frame(y = y, age = age, sex = sex, geno[, 1:200]) write.csv(model_df, "bigdata/model_demo.csv", row.names = FALSE, quote = FALSE) bm_g <- read.big.matrix( filename = "bigdata/model_demo.csv", sep = ",", header = TRUE, type = "double", backingfile = "model_demo.bin", descriptorfile = "model_demo.desc", dirpath = "bigdata/" ) fit <- bigglm( y ~ age + sex + X1 + X2 + X3, data = bm_g, family = gaussian(), chunksize = 1000 ) summary(fit)$coefficients ``` Representative coefficients on null simulated data: ``` Estimate Std. Error t value Pr(>|t|) (Intercept) 0.1234 0.2101 0.587 0.557 age -0.0021 0.0032 -0.656 0.512 sex 0.0445 0.0451 0.987 0.324 X1 0.0312 0.0478 0.653 0.514 X2 -0.0188 0.0479 -0.393 0.695 X3 0.0094 0.0477 0.197 0.844 ``` Nothing should be significant here, y is independent noise, and confirming that is the sanity check. The chunksize parameter sets how many rows bigglm reads per pass, trading speed against RAM. For SNP selection across thousands of columns, biglasso fits lasso and elastic net penalties while streaming the design matrix from disk. It requires the double type: ```r library(biglasso) bm_d <- as.big.matrix( geno, type = "double", backingfile = "geno_demo_double.bin", descriptorfile = "geno_demo_double.desc", dirpath = "bigdata/" ) X <- attach.big.matrix("bigdata/geno_demo_double.desc") fit_l <- biglasso(X, y, lambda = 0.05, penalty = "lasso") sum(coef(fit_l) != 0) ``` ``` [2] 1 ``` On null data the intercept is usually the only nonzero coefficient at this penalty strength. A char-typed matrix passed to biglasso stops with an error stating that the design matrix must be double, which is why Step 1 asks you to plan types before workflows, not during them. ### Step 9. Verify the Result and Export Three checks catch nearly everything: ```r # 1. Backing size matches the arithmetic file.size("bigdata/geno_demo.bin") == n * p * 1 # 2. Spot cells match the source all(geno[1:100, 1:50] == bm[1:100, 1:50]) # 3. Statistics agree through both paths all.equal( big.colMeans(bm, rows = 1:500, cols = 1:100), colMeans(geno[1:500, 1:100]) ) ``` ``` [2] TRUE [2] TRUE [2] TRUE ``` Export small results as ordinary R objects and plain CSV. Exporting a full matrix back to CSV with write.big.matrix works but burns time at scale, so prefer copying the backing and descriptor pair, which is a byte-for-byte backup: ```r write.big.matrix(bm_f, "bigdata/geno_maf_export.csv", sep = ",", row.names = FALSE, col.names = FALSE) ``` Archive geno_maf.bin, geno_maf.desc, and the sample key file from Step 2 together. That triple reconstructs the analysis anywhere. ## Performance Tips for Genomic Workloads 1. Orient the matrix around your dominant operation. Column-major storage makes column passes sequential. If per-sample operations dominate your workload, transpose once so samples become columns, then stream. 2. Precompute and reuse column statistics. Allele frequencies for 500,000 SNPs occupy about 4 MB as a double vector. Compute once, save to RDS, and filter from the vector rather than rescanning the matrix. 3. Keep the backing on a local solid state drive. Memory mapping over network file systems turns sequential streams into stalls, and a full pass over a 125 GB backing runs at disk read speed, which is the floor for every column statistic. 4. Match representation to information content. The two-bit genotype encodings used by biobank-scale regression tools show how much representation choice matters at this scale [1]. 5. Exploit sparsity where it exists. Rare-variant and single-cell matrices are mostly zeros, and sparse structures cut memory overhead far below any dense layout, the lesson from the CoGAPS 3 redesign [4]. 6. Factorize instead of storing. A covariance-like matrix that decomposes into a diagonal plus a low-rank part can be stored as its factors, the trick that makes correlated omics simulation extremely efficient [6]. 7. Reduce dimensions before scaling up. Principal component regression is standard in genomic prediction, and formulations such as PC-DOT, which ranks components by dot product with the response and leans on the HAT matrix to cut computation, ship with adaptable R code [7]. 8. Profile with small blocks first. Run the Step 6 loop on 100 columns, time it, then extrapolate before committing to an overnight run. Packages built for this problem advertise linear scaling with input volume for good reason [8]. ### When to Use a Dedicated Engine Instead bigmemory is a substrate for custom workflows, not a substitute for tested engines when a standard analysis exists. BESSiE implements genomic BLUP, SNP-BLUP, and the BayesA, BayesB, BayesC, and BayesR family for large-scale marker data in complex linear mixed models [2]. PRSice-2 was engineered for biobank-scale polygenic risk score analysis with memory efficiency as a design goal [9]. For phenome-wide association studies (PheWAS) on biobank data, PheTK cut runtime by 64 percent against the R PheWAS package through multithreading [10]. Even within R, population-scale problems yield to the right architecture, as LT-FGRS shows for genetic liability estimation on population-scale pedigrees at manageable computational cost [11]. Reach for bigmemory when your analysis is genuinely custom, and reach for the dedicated tool when one exists. ## Common Mistakes and Limitations 1. Extracting the whole matrix. bm[, ] and as.matrix(bm) materialize every cell in RAM. The failure is the classic allocation error, for example `Error: cannot allocate vector of size 3.7 Gb`. Work in blocks, or take a view with sub.big.matrix(bm, firstRow, lastRow), which copies nothing. 2. Choosing char for continuous dosages. Assigning 1.75 to a char cell stores 1, and the truncation is silent. After creating any matrix, spot-check a handful of cells against the source, especially fractional dosage columns. 3. Improvised missing codes. Decide the missing convention before creation, such as 3 for missing in 0, 1, 2 dosage data, write it into your codebook, and test a round trip. The double and integer types carry NA natively, and a declared missing code keeps the narrow types unambiguous whatever your bigmemory version supports. 4. Separating descriptor and backing. Renaming or moving one without the other breaks attachment, and attach.big.matrix responds with a message that the backing file cannot be found. Move the pair together, keep them in one folder, and use the backingpath argument when a relative layout changes. 5. Saving the wrong object. saveRDS on a big.matrix preserves a pointer, not data. Save describe(bm), which serializes the descriptor, and attach from it in the next session. 6. Row-wise loops on wide matrices. Iterating sample by sample across 500,000 columns re-reads scattered pages on every pass. Restructure to column blocks as in Step 6, or transpose once. 7. Parallel writes to overlapping regions. Workers reading the same backing scale cleanly. Workers writing need disjoint column or row partitions, otherwise last-writer-wins corruption lands silently in the backing. 8. Forked parallelism with a heavy parent. On macOS, forking a large R process can end the session abruptly instead of raising an error you can catch. PSOCK clusters through makeCluster avoid the pattern entirely. 9. Expecting genomics features. bigmemory stores numbers and computes basic statistics. It has no PLINK format awareness, no population genetics functions, and no plotting. Layer bigsnpr or bigstatsr for genotype-specialized operations, or dedicated engines for standard analyses. 10. Losing sample identity across sessions. Row and column names are easy to lose between attach calls. The external key file from Step 2, one ID per row in order, is the durable source of truth. ## Where to Go Next - Repeat Steps 2 through 9 on one real chromosome of dosage data, then scale to the whole genome once timings are known. - Add bigtabulate and bigalgebra if cross-tabulation or heavier matrix algebra enters your workflow. - Evaluate bigsnpr and bigstatsr if PLINK files are your native format, since they wrap the same file-backed idea with genotype-aware functions. - Benchmark against a dedicated engine such as PRSice-2 or BESSiE before promoting a custom pipeline to production. ## From Raw Genotype Files to a bigmemory Backing ![Four-step flow from VCF or PLINK files through a text export and read.big.matrix into a verified backing plus key file.](/article-images/illustrations/bigmemory-r-genomics-large-matrix-import-to-backing.webp) The demo in the main workflow starts from a matrix already in R memory, which is a luxury real pipelines rarely have. Real dosage data usually arrives as VCF, as a PLINK binary fileset, or as an exported text matrix from an imputation server. bigmemory reads delimited text directly and nothing else, so the first practical question in any project is how to get your native format into a flat, numeric, delimited file that `read.big.matrix` can stream. The pattern that works at terabyte scale has two stages: 1. Convert the native binary format to a flat delimited text file of dosages using the format's own tooling, outside R. PLINK and BCFTools both export genotype probabilities or dosages as text queries, and they do the binary parsing in compiled code far faster than anything you can write in R. 2. Stream that text file into a backing with `read.big.matrix`, keeping sample identifiers out of the matrix and in a key file. The key file deserves its own code, because it is the single habit that prevents most identity disasters downstream. If your export has samples as columns with a header row: ```r hdr <- readLines("geno_dosages.csv", n = 1) sample_ids <- strsplit(hdr, ",")[[2]][-1] # drop the leading empty or "id" cell writeLines(sample_ids, "bigdata/sample_key.txt") length(sample_ids) ``` The `-1` assumes the header begins with a label for the row-identifier column, which is the common layout. Print the first few IDs and cross-check them against the source manifest before you spend hours reading the matrix. If the export has samples as rows, read the first column of the body instead, which you can do with a connection and `readLines` in chunks so the identifier column itself never becomes a large R object. Three properties of `read.big.matrix` matter when you plan the import: - It parses the file sequentially in chunks and writes the backing as it goes, so peak RAM is governed by the chunk size, not the file size. - `type` is fixed at creation. If you declare `char` and the file contains fractional dosages, the truncation described in the main workflow happens during import, silently, across billions of cells. Inspect a sample of source lines before choosing the type. - Text parsing is the bottleneck. A pass that would take minutes from a binary backing can take hours from CSV at cohort scale, which is why the main workflow recommends binary sources when you have them. After the import completes, run the verification trio from Step 9 immediately, before any analysis: check the backing file size against the arithmetic, spot-check a block of cells against the source file, and compare a column statistic computed both ways. Import errors at scale are almost always systematic, meaning they corrupt every row or every column in the same way, and a spot check catches them in seconds while a downstream analysis can hide them for days. One further note on disk budgeting: the import temporarily needs both the text source and the backing on disk. For a whole-genome dosage export that means planning for roughly twice the text file size plus the backing size during the conversion window. Delete the intermediate text only after the verification checks pass, not before. ## A Worked Comparison: bigmemory, HDF5Array, and the bigstatsr Family The main workflow notes that Bioconductor-centric stacks tend toward DelayedArray with an HDF5 backend. Because both ecosystems appear constantly in the same project discussions, a side-by-side comparison earns its space. The differences are less about capability and more about which friction you prefer to absorb. | Dimension | bigmemory | HDF5Array / DelayedArray | bigstatsr / bigsnpr | |-|-|-|-| | Distribution | CRAN | Bioconductor | CRAN | | On-disk format | Flat binary array, uncompressed | HDF5, chunked, optional compression | Similar file-backed design with genotype-aware layer | | Sharing across processes | Native, one backing shared by many R sessions | Per-file access, no shared-memory attach | Follows bigmemory's model | | Genotype-aware operations | None, by design | None, general purpose | Yes, PLINK file input and specialized methods | | Best fit | Custom pipelines needing cross-process sharing | Bioconductor pipelines already using SummarizedExperiment | PLINK-centric genotype analysis | The compression point deserves amplification. HDF5 stores data in chunks that can be compressed, and genomic matrices compress well when they are sparse or low-entropy, so an HDF5-backed copy of a rare-variant matrix can be substantially smaller than a flat bigmemory backing. bigmemory never compresses, which is exactly why its arithmetic is so predictable: the backing size is always `nrow × ncol × bytes per cell`, and a full pass reads exactly that many bytes. Predictability is a feature when you are budgeting disk and estimating run times, and a cost when disk is the binding constraint. The sharing point matters just as much. The parallel workflow in Step 7 works because every worker attaches to one descriptor and the operating system backs all of them with the same file pages. Achieving the same effect with HDF5-backed arrays means each worker opens the file independently, which is fine for reading but loses the single-shared-copy elegance for writes. If your pipeline has a heavy parallel phase with disjoint-write patterns, bigmemory's model fits it naturally. The genotype-awareness point cuts the other way. If your native format is PLINK and your analyses are standard genotype analyses, bigsnpr and bigstatsr give you the file-backed substrate plus functions that understand allele frequencies, LD, and population structure without you writing them. The main workflow's advice stands: evaluate them first if PLINK is your native format, and reach for plain bigmemory when your analysis is genuinely custom. Whichever backend you choose, keep the interfaces between stages boring. Write per-SNP summary tables, scores, and QC flags as plain data.frames and CSV files, and keep the sample key file as the identity anchor. Backends then stay swappable: you can re-express a pipeline from bigmemory to HDF5Array by rewriting only the extraction loops, because every downstream consumer reads ordinary R objects. ## Benchmarking the Backing Before You Commit Every heavy workflow in this article reduces to the same cost model: a full pass over the backing at disk read speed, times the number of passes your algorithm needs. Measuring the cost of one pass on the demo, then extrapolating, turns overnight-run guessing into arithmetic. Time a full column pass on the demo backing: ```r system.time( afs_check <- big.colMeans(bm_f) ) ``` On a laptop this completes in well under a second, because the demo backing is only 10 MB. The useful number is throughput: ``` throughput = total bytes read / elapsed seconds ``` For the demo that is 10,000,000 bytes divided by the elapsed time. Once you have a throughput figure from your own machine, extrapolate directly: a 25 GB backing at the same throughput takes 2.5 billion divided by throughput seconds, and you can compare that against a 125 GB double-precision equivalent before anyone schedules it. Two benchmarking traps produce misleading numbers, and both are avoidable: - **Warm versus cold cache.** The operating system page cache holds recently read file pages in RAM, so a second timing pass over the same backing can be many times faster than the first because the data never touched disk. Benchmark the pass that matters, which is usually the cold one. Restart the machine or read a different large file between timings if cold behavior is what you need to know. - **Block-size sweeps that measure the wrong thing.** Very small blocks add R-level loop overhead; very large blocks can exceed your RAM budget and trigger thrashing. Sweep the middle of the range. A minimal sweep on the demo: ```r sweep_blocks <- c(100, 250, 500, 1000, 2000) timings <- sapply(sweep_blocks, function(b) { t0 <- proc.time() acc <- numeric(nrow(bm_f)) for (s in seq(1, ncol(bm_f), by = b)) { e <- min(s + b - 1, ncol(bm_f)) acc <- acc + as.numeric(bm_f[, s:e] %*% rep(1, e - s + 1)) } (proc.time() - t0)["elapsed"] }) names(timings) <- sweep_blocks timings ``` On the 10 MB demo the differences will be small and possibly noisy, which is itself the lesson: run the sweep at a scale where disk traffic dominates, for example on one real chromosome, before trusting it. The relative ordering of block sizes usually carries over to the full cohort even when absolute times do not. A companion table helps translate demo habits into cohort plans. Using the sizing formula from Step 1, at 250,000 samples with `char` storage a 2 GB block budget permits roughly 8,000 columns per block, and a 50,000-column array then needs about seven column blocks per pass. At 500,000 samples the same budget permits roughly 4,000 columns and needs about thirteen blocks. Writing these numbers down per project, next to the measured throughput, is the cheapest scheduling tool you will ever build. Finally, keep the profiling habit small and early, as the main workflow recommends: run the block loop on 100 columns first, time it, extrapolate, and only then commit the overnight run. Well-built large-data packages advertise linear scaling with input volume for good reason [8], but your pipeline is only as linear as the slowest block you wrote. ## Transposing a File-Backed Matrix the Safe Way The storage layer is column-major, so column passes stream and row passes hop. When your dominant operation is per-sample rathering than per-SNP, the main workflow's advice is to transpose once so samples become columns. This section shows the transpose itself, because doing it wrong either doubles your disk silently or crawls. The efficient pattern copies column blocks of the source into row blocks of the target: ```r n_t <- ncol(bm_f) # SNPs become rows p_t <- nrow(bm_f) # samples become columns bm_t <- big.matrix( nrow = n_t, ncol = p_t, type = "char", backingfile = "geno_t.bin", descriptorfile = "geno_t.desc", dirpath = "bigdata/" ) block <- 500 for (s in seq(1, ncol(bm_f), by = block)) { e <- min(s + block - 1, ncol(bm_f)) bm_t[s:e, ] <- t(bm_f[, s:e]) } dim(bm_t) ``` Why this ordering works: extracting columns `s:e` from the source reads a contiguous slab of the source backing, and the transpose of that slab is a compact `b × n` block that can be written into consecutive rows of the target in one assignment. Writing rows of a column-major target is strided, but writing them in blocks of 500 amortizes the stride across a large contiguous chunk of each target column. The naive alternative, looping one column at a time, turns the same job into 5,000 tiny strided writes and can be an order of magnitude slower. Verify the transpose before deleting anything: ```r all(bm_t[1:50, 1:100] == t(bm_f[1:100, 1:50])) ``` ``` [2] TRUE ``` Budget accordingly. During the transpose, both the original and the transposed backing exist on disk, so a 25 GB matrix needs about 50 GB free plus headroom. After verification, decide deliberately whether to keep both orientations: keeping both doubles disk permanently but lets every later analysis pick the orientation that streams, which at cohort scale is often worth the disk. Also remember that the transposed backing gets its own descriptor, and the sample key file now describes columns rather than rows. Document that inversion in the key file's name, for example `sample_key_transposed.txt`, or a future session will attach the wrong orientation with matching dimensions and no error. ## Computing a Genomic Relationship Matrix in Blocks The genomic relationship matrix is the canonical example of a derived object that outgrows the matrix it came from, and the main workflow quotes the arithmetic: 250,000 samples need about 500 GB at double precision before a single element is computed [3]. For cohorts up to a few thousand samples, though, a blocked in-RAM computation is entirely practical and makes an excellent teaching example, because it exercises every pattern from the main workflow at once: precomputed statistics, column-block extraction, and block accumulation. The standard GBLUP construction centers each genotype column at twice its allele frequency, forms the crossproduct of the centered matrix with itself, and scales by a function of the allele frequencies [3]. Follow the exact scaling convention of your downstream software; the blocked mechanics are identical either way. ```r afs <- big.colMeans(bm_f) / 2 n <- nrow(bm_f) G <- matrix(0, n, n) block <- 500 for (s in seq(1, ncol(bm_f), by = block)) { e <- min(s + block - 1, ncol(bm_f)) Zb <- bm_f[, s:e] # n x block, small Zb <- sweep(Zb, 2, 2 * afs[s:e]) # center each column G <- G + tcrossprod(Zb) # accumulate n x n } G <- G / (2 * sum(afs * (1 - afs))) ``` Walking the memory budget makes the pattern explicit. The extracted block `Zb` is 2,000 by 500 at double precision after extraction, a few megabytes. The accumulator `G` is 2,000 by 2000 doubles, about 32 MB, and it lives in RAM for the whole loop because every block must contribute to it. The backing is never fully resident, and peak RAM is `G` plus one block, exactly as designed. Three checks confirm the result before anyone uses it: ```r # 1. Symmetry all.equal(G[1:100, 1:100], t(G[1:100, 1:100])) # 2. Diagonal centered near 1 for the standard scaling summary(diag(G)) # 3. Off-diagonals near 0 for unrelated simulated samples off <- G[upper.tri(G[1:200, 1:200])] summary(off) ``` On simulated unrelated samples the diagonal should cluster near 1 and the off-diagonal pairs near 0, with small positive values from sampling noise. Real cohorts break this pattern in informative ways: family structure, population stratification, and duplicated samples all show up as diagonal or off-diagonal outliers, which is why a quick summary of `diag(G)` and a sample of off-diagonal pairs is a QC step in its own right. Now the honest limitation. The accumulator `G` scales with the square of sample size, so the pattern that works at 2,000 samples fails hard at 250,000, where `G` alone demands the 500 GB quoted earlier. At that scale the choices are a low-rank approximation, since principal-component formulations that lean on the HAT matrix reduce exactly this computational burden [7], factorized storage of a diagonal-plus-low-rank structure in the style that makes correlated omics simulation efficient [6], or a dedicated engine built for genomic prediction at scale [2]. The blocked crossproduct remains valuable below that threshold and, just as importantly, it is the template for any statistic that reduces a tall-skinny matrix to a small square object. ## Missing Data Strategies That Stay in Bounds Missingness handling is where the narrow element types earn and demand care. The main workflow fixes the convention of 3 for missing in 0, 1, 2 dosage data and shows how `mwhich` locates such cells. This section covers what to do with them, in bounded memory, and what each storage type permits. First, the type facts, because they decide your options: - `double` and `integer` carry `NA` natively. A cell can be genuinely missing and `mwhich` and column statistics handle it. - `char` and `short` carry no `NA`. Every cell is an integer, so missingness must be encoded as a reserved value such as 3, and you must remember that the value is a code, not a measurement. Second, the storage consequence of imputation. Mean imputation sets a missing dosage to twice the allele frequency, for example 0.75 for an allele frequency of 0.375. That value cannot live in a `char` cell: it truncates to 0, corrupting the matrix in exactly the silent way the main workflow warns about. Two honest options exist: 1. **Round imputed values to valid dosages.** Impute 0, 1, or 2 by nearest dosage. This keeps `char` storage but distorts the imputation, and the distortion interacts badly with allele-frequency-based centering downstream. 2. **Scale and widen.** Store dosages multiplied by a fixed factor in a wider type. With `short`, which holds integers up to about 32,767, a scale factor of 1,000 stores dosages from 0 to 2,000 with three decimal places of headroom, and the article's type table already flags this scaled-dosage use. The scaled option, in blocks: ```r bm_imp <- big.matrix( nrow = nrow(bm_f), ncol = ncol(bm_f), type = "short", backingfile = "geno_imp.bin", descriptorfile = "geno_imp.desc", dirpath = "bigdata/" ) block <- 500 for (s in seq(1, ncol(bm_f), by = block)) { e <- min(s + block - 1, ncol(bm_f)) B <- bm_f[, s:e] is_miss <- B == 3 for (j in seq_len(ncol(B))) { B[is_miss[, j], j] <- 2 * afs[s + j - 1] } bm_imp[, s:e] <- round(B * 1000) # scaled encoding, decode by /1000 } ``` The inner loop over columns of the block is short, since the block is only 500 columns wide, and the outer loop keeps peak RAM at the block size. The scale factor must be written into the project codebook next to the missing-code convention, because a scaled `short` backing is unreadable six months later without it. A reader who decodes by dividing by 1,000 when the writer used 10,000 produces a uniformly wrong matrix with plausible-looking statistics, one of the most dangerous failure modes in file-backed work because every downstream check can pass. There is also a third option worth naming: do not imitate at all. Many column-wise analyses tolerate missing codes if you handle them explicitly. In the relationship-matrix computation above, a missing cell that you replace by the column mean before centering contributes exactly zero to the crossproduct, which is equivalent to mean imputation without writing a second backing. Choosing per analysis, rather than baking one imputation into the archive copy, keeps your options open and your backing honest. ## Parallel Patterns That Scale Past Column Means Step 7 shows the canonical read-only pattern, workers attaching to one descriptor and consuming disjoint column ranges. Three patterns extend that idea to the analyses genomics actually runs, and each has a verification habit attached. **Pattern 1: partitioned model fits.** When a genome-wide scan decomposes by chromosome, give each worker a chromosome's column range and let it fit whatever per-marker or per-region model the analysis specifies. Because chromosomes are contiguous column ranges in a properly ordered backing, each worker reads a contiguous slab and disk behavior stays sequential: ```r chrom_bounds <- list(c(1, 1500), c(1501, 3200), c(3201, 4700), c(4701, 5000)) cl <- makeCluster(length(chrom_bounds)) registerDoParallel(cl) results <- foreach( bounds = chrom_bounds, .packages = "bigmemory", .combine = "rbind" ) %dopar% { bm_w <- attach.big.matrix(desc) chunk <- bounds[2]:bounds[10] afs_w <- big.colMeans(bm_w, cols = chunk) data.frame(snp_index = chunk, af = afs_w / 2) } stopCluster(cl) nrow(results) ``` The `.combine` stitching depends on the column ranges being disjoint and sorted, which the bounds list guarantees. Assert that invariant in code rather than trusting the list: overlapping ranges produce silently wrong parallel results, the same class of failure as overlapping writes. **Pattern 2: row-partitioned cross-validation.** Views created with `sub.big.matrix` copy nothing, so k-fold cross-validation can hand each fold a view without extracting data. Workers attach the descriptor, build their fold's row view, and fit or score independently: ```r folds <- cut(seq_len(nrow(bm_f)), breaks = 5, labels = FALSE) fold_results <- foreach( k = 1:5, .packages = "bigmemory", .combine = "c" ) %dopar% { bm_w <- attach.big.matrix(desc) rows_k <- which(folds == k) view <- sub.big.matrix(bm_w, firstRow = min(rows_k), lastRow = max(rows_k)) mean(big.colMeans(view)) # replace with the real fold statistic } ``` This simplified example computes a trivial statistic; the structural point is that the view is a pointer, not a copy, so fold extraction costs nothing and the row assignment stays honest because `cut` produces contiguous blocks on the demo. On real data with shuffled sample order, sort the fold assignment first or accept strided reads. For modeling work inside folds, remember the type constraints from Step 8: penalized fitters that require `double` reject `char` backings at fit time. **Pattern 3: disjoint-write partitioning.** The main workflow warns that overlapping parallel writes corrupt the backing silently. The constructive version of that warning is to design imputation and filtering jobs so that each worker owns a column range exclusively and writes only to its own columns of a preallocated target backing. The transposition and imputation loops above are already shaped this way, so parallelizing them is a mechanical change: give each worker a start index, have it process and write its range, and keep the column-range table as the audit record of who wrote what. A word on architecture choice, because it frames when these patterns beat dedicated tools. Multithreading inside a well-engineered application delivered a 64 percent runtime reduction for phenome-wide association work relative to an R package [10], and population-scale pedigree analyses have been made tractable in R through architecture rather than hardware alone [11]. The lesson for bigmemory users is that the parallel patterns above buy real speed on custom analyses, but a mature dedicated engine with threading baked in usually wins on standard analyses, as the main workflow's guidance already concludes. Every parallel rewrite in this article ends with the same habit: compare the parallel result against a serial pass with `all.equal` on a bounded subset. That check costs seconds and catches the partition arithmetic errors that otherwise reach production silently. ## Sampling, Subsetting, and Permutation From Disk Large backings change the economics of sampling. Extracting a random subset of rows or columns is cheap precisely because it touches few pages, and several routine tasks reduce to exactly that. **Random row samples for exploration.** Sampling a few hundred rows and pulling them fully into RAM as an ordinary matrix opens the door to every standard R tool, plots included: ```r set.seed(11) row_sample <- sample(nrow(bm_f), 200) sub_mat <- bm_f[row_sample, 1:200] # 200 x 200, trivially small dim(sub_mat) summary(colMeans(sub_mat)) ``` Two disciplines keep this honest. Record the sampled indices to an RDS file so any figure made from the subset is reproducible, and remember that row sampling from a column-major backing reads strided pages, so a 200-row sample of 5,000 columns touches the whole file lightly rather than reading a contiguous slab. That is fine at 10 MB and worth thinking about at 25 GB. **Column bootstraps for uncertainty.** Bootstrap resampling of markers, with replacement, is a column-index operation plus a block pass: ```r set.seed(21) boot_idx <- sample(ncol(bm_f), size = ncol(bm_f), replace = TRUE) boot_afs <- big.colMeans(bm_f, cols = boot_idx) / 2 quantile(boot_afs, c(0.025, 0.975)) ``` Because duplicate indices simply re-read the same pages, which the page cache serves quickly, this costs little more than a single pass. The same trick produces bootstrap distributions for any column-wise statistic you can express in a block loop, and it pairs naturally with the precomputed-statistics habit from the performance tips: resample the saved allele-frequency vector for free, and only resample the backing when the statistic requires the raw genotypes. **Permutation structure for significance.** Permutation testing reshuffles the response and repeats a statistic, and the disk-backed version streams the matrix once per permutation. The cost model from the benchmarking section applies directly: a permutation analysis at B replicates costs roughly B passes, so estimate the pass time first and let it decide B, rather than picking B from convention and discovering the runtime later. For genome-wide significance testing at biobank scale, dedicated tools remain the pragmatic choice since they were engineered with exactly this cost model in mind [9], and custom permutation loops are best reserved for bespoke statistics that no engine computes. One caution ties this section back to study design. Sampling rows before model fitting changes the population being modeled, and marker-level sampling before a scan changes the search space. Use sampling for exploration, uncertainty, and testing, and keep the full-cohort steps full-cohort, with the block patterns doing the bounded-memory work. ## Troubleshooting: Symptoms, Causes, and Fixes ![Decision path from allocation errors, attachment failures, wrong numbers and parallel mismatches to the minimal fix.](/article-images/illustrations/bigmemory-r-genomics-large-matrix-troubleshoot-triage.webp) Most bigmemory failures cluster into a small set of recognizable shapes. Matching symptom to cause quickly is the difference between a five-minute fix and a lost afternoon, so this section pairs each common symptom with its usual cause and the minimal fix. | Symptom | Most likely cause | Fix | |-|-|-| | `Error: cannot allocate vector of size X Gb` | A whole-matrix extraction such as `bm[, ]`, `as.matrix(bm)`, or an unblocked assignment hid in the code | Find the unguarded extraction, replace with a block loop or a `sub.big.matrix` view | | Attachment fails with a message that the backing file cannot be found | Descriptor and backing moved or renamed separately | Move the pair together into one folder, or supply the correct `backingpath` | | Numbers look wrong only in some columns | `char` truncation of fractional values, or a scaled-encoding factor mismatch | Spot-check those columns against the source, re-import with a wider type if truncation is confirmed | | `biglasso` stops immediately complaining about matrix type | A `char` backing passed where `double` is required | Create or copy the matrix as `double` for the modeling stage, as Step 8 shows | | Parallel result differs from the serial result | Overlapping column ranges between workers, or a missing `attach.big.matrix` inside the worker | Print the range table, assert disjointness, confirm every worker attaches inside `%dopar%` | | Values vanish or change after a parallel run | Two workers wrote overlapping regions, last writer won | Redesign to disjoint-write partitions, restore from the verified backup | | Session ends abruptly with no R error under parallel load | Forked parallelism from a memory-heavy parent on macOS | Switch to `makeCluster` PSOCK workers | | Second run of the same code is mysteriously much faster | Page cache serving previously read pages | Not an error, but do not benchmark on warm cache when cold behavior matters | | Row or column identity scrambled after reattachment | Identifier handling drifted between sessions | Rebuild from the key file, one ID per row in file order, and assert dimensions before use | | Backing file exists but seems smaller than expected | Creation-time truncation, often fractional dosages into `char` | Recompute `nrow × ncol × bytes per cell` and compare against `file.size` before trusting anything | Three of these deserve extra emphasis because they do not fail loudly. Truncation fails silently in every cell of every affected column. Overlapping parallel writes fail silently in whichever regions raced. And warm-cache speedups fail silently in your planning arithmetic, producing run-time estimates you can never reproduce on a cold cluster node. The main workflow's verification trio, backing size arithmetic, spot cells, and dual-path statistics, catches the first two; the benchmarking section's cold-cache discipline catches the third. Run all three checks as a scripted preflight rather than as memory. When a failure does reach you as a corrupted backing, the recovery path is the archive habit from Step 9: the backing, descriptor, and key file triple. Backings are immutable in well-designed pipelines, since analyses read them and write their outputs elsewhere, which means the backup from the last verified state is always sufficient. Pipelines that write in place to production backings lose this safety net, and that design choice is usually visible in the first incident report. ## Misconceptions Worth Retiring Several durable misconceptions about bigmemory circulate in genomics training materials and code review comments. Correcting them early prevents misdesigned pipelines. **Misconception 1: memory mapping means RAM-speed access.** The operating system maps file pages into address space, but the first touch of each page still comes from disk, and eviction under memory pressure sends performance toward disk speed. Bigmemory removes the size ceiling, not the storage hierarchy. Design for sequential column passes, as the whole main workflow does, and the disk speed you get is sequential-read speed, which is the good kind. **Misconception 2: the backing is a database.** It is a flat array of numbers with fixed dimensions and no index, no query language, and no compression. Anything resembling a query, such as finding rows where several conditions hold, is a scan you write yourself with tools like `mwhich`, and its cost is a pass. **Misconception 3: `char` means text.** In base R, `char` might suggest character strings. In bigmemory, `char` is a one-byte signed integer, which is why dosages fit and why sample names do not. The article's type table exists largely to prevent this confusion at creation time, when it is cheap to fix. **Misconception 4: attaching copies the data.** `attach.big.matrix` builds a small object holding a pointer and metadata. A hundred sessions can attach to one backing and the disk holds one copy. This is the foundation of the parallel patterns, and forgetting it leads people to avoid parallel workflows that are in fact the cheapest thing in the toolbox. **Misconception 5: bigmemory handles sparse data.** It stores dense arrays, period. For mostly-zero matrices the correct move is a sparse representation, and the published evidence is emphatic about how much representation matters: two-bit genotype encoding makes biobank-scale regression feasible where eight-byte dense storage would not [1], and sparse structures carried a Bayesian factorization workflow to roughly a thousandfold more cells [4]. Dense file backing is a choice for dense data, not a default for all data. **Misconception 6: file-backed means crash-safe.** The backing reflects writes as they happen, but an interrupted job can leave a half-written target backing exactly as an interrupted write to any file can. Crash safety comes from the pipeline discipline of writing to new files, verifying, and only then replacing archived versions, not from the storage layer. **Misconception 7: once it is on disk, statistics are free.** Every column statistic is a full pass at disk speed. Precomputing the 4 MB allele-frequency vector once, as the performance tips advise, is the difference between one pass per project and one pass per question. ## Backing File Hygiene and Reproducible Pipelines File-backed matrices change what reproducibility means, because the primary data artifact is now a binary pair plus a key file rather than a script's intermediate object. A small amount of hygiene makes the whole arrangement durable. Record integrity at creation, immediately after the Step 9 checks pass: ```r tools::md5sum(c("bigdata/geno_demo.bin", "bigdata/geno_demo.desc")) ``` Store the hash values in a plain text manifest next to the files. Any later pipeline stage can recompute them and refuse to proceed on mismatch, which converts silent corruption from a downstream mystery into an immediate, local error. For backings that a pipeline regenerates, compute the hash at the end of the generating script and have consuming scripts verify it as their first action. Write the provenance record the same way. A backing without provenance is an unlabeled freezer sample; a few lines of structured text prevent that: ```r provenance <- c( "backing = geno_demo.bin", "descriptor = geno_demo.desc", "key_file = sample_key.txt", "nrow = 2000", "ncol = 5000", "type = char", "coding = 0/1/2 dosage, 3 = missing", "scale_factor = 1 (unscaled)", "source = simulated, set.seed(42)", "created = 2026-11-11", "verified_against_source = TRUE" ) writeLines(provenance, "bigdata/geno_demo_README.txt") ``` Adapt the fields to the project, but keep three nonnegotiable: the coding convention, the scale factor if any, and the verified-against-source flag. Every troubleshooting entry in the previous section becomes easier, and several become unnecessary, when those three facts travel with the file. Capture the software context with the same discipline: ```r sessionInfo() # save its output to bigdata/sessionInfo.txt in the same script that # creates the backing ``` Package versions affect file-backed behavior across the ecosystem, and a backing created under one combination may attach under another with different edge-case behavior. The `sessionInfo` snapshot costs nothing and answers the question every incident investigation asks first. Finally, treat version control sensibly. Backings are large binaries and do not belong in a code repository; the descriptor, key file, provenance record, and hash manifest are tiny and absolutely do belong there. A pipeline repository that contains every descriptor and key file, plus the scripts that regenerate backings from sources, reconstructs the entire analysis state from any checkout, with the backings themselves restorable from archives. ## Pre-Flight Checklist for an Overnight Run Before committing hours of compute to a large file-backed job, walk this list. Each item traces to a section above, and each has caught real failures. 1. **Sizing arithmetic written down.** Backing bytes computed from `nrow × ncol × bytes per cell`, compared against free disk, with room for both orientations if a transpose is planned. 2. **Type chosen against evidence.** Source values inspected for fractions; missing-code convention documented; scale factor documented if `short` scaling is in play. 3. **Cold-pass throughput measured.** One real pass timed with a cold page cache, run time for the full job extrapolated from the measured throughput. 4. **Block size chosen from a sweep, not a guess.** Sweep run at a meaningful scale, block size recorded in the project notes next to the RAM budget that produced it. 5. **Extraction audit clean.** No unguarded `bm[, ]`, `as.matrix(bm)`, or unblocked column assignment anywhere in the job script. 6. **Parallel partition table printed and asserted disjoint.** Worker ranges verified nonoverlapping, every worker attaches inside the parallel block, serial-vs-parallel `all.equal` check scripted. 7. **Write targets disjoint or serial.** Any parallel write confirmed to touch exclusive column or row ranges, or the stage downgraded to serial. 8. **Verification trio scripted at the end of the job.** Backing size arithmetic, spot cells against source, and dual-path statistics, all automated, all failing loudly. 9. **Manifest complete.** Hashes computed, provenance written, key file present, `sessionInfo` captured, descriptor and key file committed to the repository. 10. **Rollback defined.** The previous verified backing pair retained until the new one passes checks, and the job writes new files rather than overwriting archived ones. The checklist costs minutes to complete and encodes every failure mode this guide has documented, from silent truncation through last-writer-wins corruption to misplanned run times. As cohort sizes continue climbing and matrices that once demanded clusters approach laptop scale through good representation choices [1], the workflows stay the same even as the numbers change, and the habit that transfers best is exactly this one: verify the storage layer before trusting anything computed on it, whether the analysis is a custom blocked crossproduct, a penalized fit streamed from disk, or a genomic prediction that could equally have gone to a dedicated engine built for the purpose [2]. ## Frequently Asked Questions ### How large can a bigmemory matrix get? ## Sources 1. [Fast numerical optimization for genome sequencing data in population biobanks.](https://pubmed.ncbi.nlm.nih.gov/34146108/) 2. [BESSiE: a software for linear model BLUP and Bayesian MCMC analysis of large-scale genomic data.](https://pubmed.ncbi.nlm.nih.gov/27590176/) 3. [Sharing of either phenotypes or genetic variants can increase the accuracy of genomic prediction of feed efficiency.](https://pubmed.ncbi.nlm.nih.gov/36068488/) 4. [CoGAPS 3: Bayesian non-negative matrix factorization for single-cell analysis with asynchronous updates and sparse data structures.](https://pubmed.ncbi.nlm.nih.gov/33054706/) 5. [GWAtoolbox: an R package for fast quality control and handling of genome-wide association studies meta-analysis data.](https://pubmed.ncbi.nlm.nih.gov/22155946/) 6. [Generating correlated data for omics simulation.](https://pubmed.ncbi.nlm.nih.gov/40911660/) 7. [PC-DOT: Improving genomic prediction ability of principal component regression by DOT product.](https://pubmed.ncbi.nlm.nih.gov/36168679/) 8. [minSNPs: an R package for the derivation of resolution-optimised SNP sets from microbial genomic data.](https://pubmed.ncbi.nlm.nih.gov/37250706/) 9. [PRSice-2: Polygenic Risk Score software for biobank-scale data.](https://pubmed.ncbi.nlm.nih.gov/31307061/) 10. [PheWAS analysis on large-scale biobank data with PheTK.](https://pubmed.ncbi.nlm.nih.gov/39657951/) 11. [LT-FGRS: a unifying R-package for the estimation of family-based genetic liabilities at population-scale.](https://pubmed.ncbi.nlm.nih.gov/42367923/)