# Read and write BED files with genio In genio: Genetics Input/Output Functions

knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )  ## copied from examples from the "simmer" R package ## after: https://www.enchufa2.es/archives/suggests-and-vignettes.html ## by Iñaki Úcar required <- c("BEDMatrix", "snpStats", "pryr") if (!all(sapply(required, requireNamespace, quietly = TRUE))) knitr::opts_chunk$set(eval = FALSE)


# Introduction

The genio (GenIO = Genetics I/O) package aims to facilitate reading and writing genetics data. The focus of this vignette is processing Plink BED/BIM/FAM files.

There are some limited alternatives for reading and/or writing BED files in R, which are slower and harder to use, which motivated me to write this package. Here we make direct comparisons to those packages, to illustrate the advantages of genio.

# Generate a large matrix of random genotypes

Let's begin by creating a large genotype matrix with completely random data, and with missing values so this becomes non-trivial to read and write.

# Data dimensions.
# Choose non-multiples of 4 to test edge cases of BED parsers.
# Number of loci.
m_loci <- 10001
# Number of individuals
n_ind <- 1001
# Overall allele frequency
# (we don't do any actual genetics in this vignette,
# so just set to a reasonable value)
p <- 0.5
# Missingness rate
miss <- 0.1

# Total number of genotypes
n_data <- n_ind * m_loci
# Draw random genotypes from Binomial
X <- rbinom( n_data, 2, p)
X[ sample(n_data, n_data * miss) ] <- NA
# Turn into matrix
X <- matrix(X, nrow = m_loci, ncol = n_ind)

# Inspect the first 10 individuals at the first 10 loci
X[1:10, 1:10]


To create annotation tables that look a bit more interesting than the defaults, let us create some slightly more realistic values. Here we use our first genio function!

library(genio)


First we create and edit the locus annotations table.

# We have to specify the number of loci
bim <- make_bim( n = m_loci )

# Inspect the default values
bim

# Let's add the "chr" prefix to the chromosome values,
# so we recognize them when we see them later.
bim$chr <- paste0('chr', bim$chr)
# Make SNP IDs look like "rs" IDs
bim$id <- paste0('rs', bim$id)
# Make positions 1000 bigger
bim$pos <- bim$pos * 1000
# Select randomly between Cs and Gs for the reference alleles
bim$ref <- sample(c('C', 'G'), m_loci, replace = TRUE) # Select randomly between As and Ts for the alternative alleles bim$alt <- sample(c('A', 'T'), m_loci, replace = TRUE)

# Inspect the table with our changes
bim


Now we similarly create and edit the individual annotations table.

# Specify the number of individuals
fam <- make_fam( n = n_ind )

# Inspect the default values
fam

# Add prefixes to families and IDs to recognize them later
fam$fam <- paste0('fam', fam$fam)
fam$id <- paste0('id', fam$id)
# Sex values are usually 1 and 2
fam$sex <- sample(1:2, n_ind, replace = TRUE) # Let's make phenotypes continuous. # Draw independently from Standard Normal. fam$pheno <- rnorm(n_ind)
# Let's leave maternal and paternal IDs as missing

# Inspect again
fam


Lastly, let's copy the locus and individual IDs as row and column names of the genotype matrix, respectively. Although this step is not required, it is encouraged for consistency (if present, values are checked when writing files).

# Add column and row names from bim and fam tables we just created.
rownames(X) <- bim$id colnames(X) <- fam$id
# Inspect again the first 10 individuals and loci
X[1:10, 1:10]


# Write genotypes to Plink files

Let's write this random data to a file. This mode is intended for simulated data, generating dummy BIM and FAM files in the process and writing out all three.

# Will delete at the end of the vignette

# Write genotypes, along with the BIM and FAM files we created.
# Omiting them would result in writing the original dummy version of these tables,
# before we edited them.
time_write_genio <- system.time(
)
time_write_genio


Here we demonstrate how easy it is to read the data back. To compare to other packages, we shall time loading all this data.

# Read the data back in memory.
# Time this step
)

# Inspect the data we just read back

# The same random genotypes (first 10 individuals and loci, now with row and column names):
data_genio$X[1:10, 1:10] # The locus annotations data_genio$bim

# The individual annotations
data_genio$fam # Quickly test that the inputs and outputs are identical. # Genotypes have NAs, so we have to compare this way. stopifnot( all( X == data_genio$X, na.rm = TRUE) )
stopifnot( bim == data_genio$bim ) # FAM has mixed content (chars, ints, and doubles). # First 5 columns should be exact match: stopifnot( fam[,1:5] == data_genio$fam[,1:5] )
# Exact equality may fail for pheno due to precisions, so test this way instead:
stopifnot( max(abs(fam$pheno - data_genio$fam$pheno)) < 1e-4 )  # Memory estimation One drawback of genio and other related approaches is high memory consumption (see more on that at the end of this vignette). It is entirely possible that an attempt to parse a genotype matrix into memory will fail with an "out of memory" error message. Let's estimate memory usage here. Genotypes are stored as integers by genio, which in R on a modern machine (64 bit architecture) consumes 4 bytes! So an$n \times m$matrix takes up at least$4 n m$bytes (an R matrix contains an additional constant overhead, which does not depend on$n,m$and is relatively small for large matrices). To get an idea of how much this is, let's assume a standard genotyping array, where$m$is about half a million. In this scenario, we get a simple rule of thumb that every 1000 individuals consume a little less than 2G: # Constants bytes_per_genotype <- 4 bytes_per_gb <- 1024 ^ 3 # Example data dimensions num_ind <- 1000 num_loci <- 500000 # Gigabytes per 1000 individuals for a typical genotyping array bytes_per_genotype * num_ind * num_loci / bytes_per_gb  This is the amount of free memory required just to load the matrix. Several common matrix operations for genotypes consume even more memory, so more free memory will be needed to accomplish anything. # Comparison to other packages There are several R packages that provide BED readers, and to a lesser extent some of the other functionality of genio. Here we compare directly to BEDMatrix and snpStats. Each of these packages has different goals so they are optimized for their use cases. The genio parser is optimized to read the entire genotype data into a native R matrix, so that it is easy to inspect and manipulate for R beginners. Here we demonstrate that genio is not only the fastest at this task, but also the easiest to obtain in terms of requiring the least amount of coding. ## BEDMatrix The BEDMatrix package allows access to genotypes from a BED file without loading the entire matrix in memory. This is very convenient for large datasets, as long as only small portions of the data are pulled at once. However, there are some disadvantages: • The BEDMatrix return object is not a regular R matrix, which confuses some users. • The BEDMatrix package provides no way to access the full annotation tables (BIM and FAM files). • From BIM, only the locus ID and reference allele are returned, merged into a vector of strings and placed as column names. • From FAM, only the family and individual IDs are returned, merged into a vector of strings and placed as row names. • BEDMatrix does not have any write functions. Here is an example of that usage and how the data compares. Note in particular that the BEDMatrix representation of the genotype matrix is transposed compared to the genio matrix: library(BEDMatrix) # Time it too. # Although the BIM and FAM tables are not returned, # they are partially parsed and kept in memory, # which can take time for extremely large files time_read_bedmatrix_1 <- system.time( X_BEDMatrix <- BEDMatrix(file_plink) ) # Inspect the first 10 loci and individuals as usual. # Note it is transposed compared to our X. # Also note the column and row names are different from genio's. X_BEDMatrix[1:10, 1:10]  Therefore, for very large datasets, if it suffices to access the genotype data in small portions and the user is willing to deal with the limitations of the object, and no detailed annotation table information is required, then BEDMatrix is a better solution than the read_bed function provided in the genio package. To compare the two genotype reading functions on an equal footing, let us assume that the entire genotype matrix is required in memory and stored in a regular R matrix. # This turns it into a regular R matrix. # Since most of the reading is actually happening now, # we time this step now. time_read_bedmatrix_2 <- system.time( X_BEDMatrix_Rmat <- as.matrix(X_BEDMatrix) ) time_read_bedmatrix_2 # Now we can test that the BEDMatrix agrees with the original matrix we simulated. # Need to transpose first. stopifnot( all( X == t(X_BEDMatrix_Rmat), na.rm = TRUE) )  ## snpStats The snpStats package is the most fully-featured alternative to genio. One of its advantages is its memory efficiency, encoding the entire data in less memory than a native R matrix. However, this memory efficiency comes at the cost of making the data harder to access, especially for inexperienced R users. Like genio, and in contrast to the other packages, snpStats also reads the BIM and FAM tables fully, and provides a BED writer, but it is considerably harder to use. First we illustrate data parsing. The annotation tables are similar but column names are different, and certain missing values (zero in text files) are converted to NA instead. library(snpStats) # Read data, time it. time_read_snpStats_1 <- system.time( data_snpStats <- read.plink(file_plink) ) time_read_snpStats_1 # Inspect the data # Genotypes. # Note data is not visualized this way. # This matrix is also transposed compared to the genio matrix. data_snpStats$genotypes

# Locus annotations
head( data_snpStats$map ) # Individual annotations head (data_snpStats$fam )


As for BEDMatrix, assuming we ultimately desire to convert the entire data into a regular R matrix, an extra step is required:

# Transpose, then convert to a regular R matrix.
# Let's time this step too.
X_snpStats <- as( t(data_snpStats$genotypes), 'numeric') ) time_read_snpStats_2 # Now we can visualize the matrix. # First 10 loci and individuals, as before. # Note that, compared to (genio, BEDMatrix), alleles are encoded in reverse, # so 0s and 2s are flipped in this matrix. X_snpStats[1:10, 1:10] # Again verify that the matrices are identical. # (Here 2-X flips back 0s and 2s) stopifnot( all( X == 2 - X_snpStats, na.rm = TRUE) )  snpStats is the only package other than genio to provide a BED writer. Here we demonstrate how hard it is to use it to write our data. # Let's write this to another file file_plink_copy <- tempfile('vignette-random-data-copy') # Copy objects to not change originals X_snpStats <- X bim_snpStats <- as.data.frame(bim) # to use rownames fam_snpStats <- as.data.frame(fam) # ditto # All data requires matching row and/or column names. # These first two were already done above. #rownames(X_snpStats) <- bim$id
#colnames(X_snpStats) <- fam$id # Row names here are redundant but required. rownames(bim_snpStats) <- bim$id

# Create barplot

barplot( sizes, names.arg = names_sizes, main = 'Native genotype object sizes', xlab = 'packages', ylab = 'memory (bytes)' )

The BEDMatrix package is ideal for very large files, since data is loaded into memory only as needed.
So the main object does not actually hold any genotypes in memory.
It is up to the user to decide how much data to access at once to trade-off speed and memory consumption.

snpStats manages memory better than genio, but strictly by a factor of 4.
Each genotype is encoded as an integer in R (in general) and genio (in particular), which uses 4 bytes (this is actually platform dependent).
On the other hand, snpStats uses strictly one byte per genotype.

Overall, there is a narrow window in data sizes in which a genotype matrix is too large for a native R matrix but small enough for a snpStats object.
That, and the much more complex interface of snpStats, makes it not very worthwhile in my opinion, unless there is need for functions provided only by that package.
I prefer to use BEDMatrix for large datasets.

# Cleanup

This handy function removes the three BED/BIM/FAM files we generated at the beginning of this vignette.

{R}

sessionInfo()