Genomic data analysis with large matrices in R

Example of expression data

Let’s use the NCBI GEO archive to get some example expression data (Barrett et al. 2012).

# http://bioconductor.org/packages/release/bioc/html/GEOquery.html
suppressMessages(library(GEOquery))

# All data
gds <- getGEO("GDS3929")

## File stored at:

## /tmp/RtmpOTdb2V/GDS3929.soft.gz

str(gds, max.level = 2)

## Formal class 'GDS' [package "GEOquery"] with 3 slots
##   ..@ gpl      :Formal class 'GPL' [package "GEOquery"] with 2 slots
##   ..@ dataTable:Formal class 'GEODataTable' [package "GEOquery"] with 2 slots
##   ..@ header   :List of 23

format(object.size(gds), units = "MB")

## [1] "37.1 Mb"

df <- gds@dataTable@table
format(object.size(df), units = "MB")

## [1] "37 Mb"

# Part of the data that can be encoded as a matrix
mat <- t(df[-(1:2)])
dim(mat)

## [1]   183 24526

format(object.size(mat), units = "MB")

## [1] "35.6 Mb"

In this data, 96% of the information can be encoded as a matrix. For larger data, this could increase to more than 99%. Let us store this as a Filebacked Big Matrix (FBM).

# devtools::install_github("privefl/bigstatsr")
library(bigstatsr)
fbm <- as_FBM(mat)
# Same data?
all.equal(fbm[], mat, check.attributes = FALSE)

## [1] TRUE

# Size in memory
format(object.size(fbm), units = "MB")

## [1] "0 Mb"

# File where the data is backed
fbm$backingfile

## [1] "C:\\Users\\Florian\\AppData\\Local\\Temp\\Rtmps7eEQV\\file1bbc9c85c43.bk"

So, now the matrix is stored on disk and basically takes no memory.

For omic data, computing PCA (or SVD) is sometimes useful to assess the structure in the data. Let’s do this.

system.time(
  svd <- svd(scale(mat), nu = 10, nv = 10)
)

##    user  system elapsed 
##   2.265   0.042   2.316

plot(svd$u)

system.time(
  svd2 <- big_SVD(fbm, big_scale(), k = 10)
) # see big_randomSVD() for very large data

##    user  system elapsed 
##    0.74    0.12    0.88

plot(svd2)

plot(svd2, type = "scores")

The goal of memory-mapping is to not have all the data in memory. This force programmers to think more before implementing algorithms. This can lead to implementations that can be much faster than current R implementations. Moreover, data stored on disk is shared between processes so that parallelism in algorithms can be implemented more easily and efficiently, for example removing the need to copy the data to each core.

Example of SNP data

Let’s use some data from Gad Abraham’s GitHub repository that have been filtered and quality-controled.

tmpfile <- tempfile()
base <- paste0(
  "https://github.com/gabraham/flashpca/raw/master/HapMap3/",
  "1kg.ref.phase1_release_v3.20101123_thinned_autosomal_overlap")
exts <- c(".bed", ".bim", ".fam")
purrr::walk2(paste0(base, exts), paste0(tmpfile, exts), ~download.file(.x, .y))

First, we have to transform these data in PLINK format to an FBM. In package {bigsnpr}, function snp_readBed() is provided to do just that. It is implemented with Rcpp, which makes it very fast. Package {Rcpp} makes easy to use C++ code in R and with R objects.

# https://github.com/privefl/bigsnpr
library(bigsnpr)
rdsfile <- snp_readBed(paste0(tmpfile, ".bed"))

## Warning in snp_readBed(paste0(tmpfile, ".bed")): EOF of bedfile has not been reached.

bigsnp <- snp_attach(rdsfile)
(G <- bigsnp$genotypes)

## A Filebacked Big Matrix of type 'code 256' with 1092 rows and 14079 columns.

object.size(G)

## 688 bytes

G[which(is.na(G[]))] <- 0  # quick and dirty imputation
G.matrix <- G[]
format(object.size(G.matrix), units = "MB")

## [1] "117.3 Mb"

For large SNP data, this format can save us a lot of memory.

system.time(
  svd <- svd(scale(G.matrix), nu = 10, nv = 10)
)

##    user  system elapsed 
##   9.154   0.741   9.895

plot(svd$u)

system.time(
  svd2 <- big_SVD(G, big_scale(), k = 10)
)

## (2)

##    user  system elapsed 
##   0.941   0.045   0.986

plot(svd2)

plot(svd2, type = "scores")