R Interface to CoreArray Genomic Data Structure (GDS) Files
Xiuwen Zheng
Jul 5, 2015
Introduction
The package provides a high-level R interface to CoreArray Genomic Data Structure (GDS) data files, which are portable across platforms and include hierarchical structure to store multiple scalable array-oriented data sets with metadata information. It is suited for large-scale datasets, especially for data which are much larger than the available random-access memory. The package gdsfmt offers the efficient operations specifically designed for integers with less than 8 bits, since a single genetic/genomic variant, like single-nucleotide polymorphism (SNP), usually occupies fewer bits than a byte. Data compression and decompression are also supported with relatively efficient random access.
Installation of the package gdsfmt
To install the package gdsfmt, you need a current version (>=2.14.0) of R. After installing R you can run the following commands from the R command shell to install the package gdsfmt.
Install the package from Bioconductor repository:
source("http://bioconductor.org/biocLite.R")
biocLite("gdsfmt")Install the development version from Github:
library("devtools")
install_github("zhengxwen/gdsfmt")The install_github() approach requires that you build from source, i.e. make and compilers must be installed on your system – see the R FAQ for your operating system; you may also need to install dependencies manually.
High-level R functions
Creating a GDS file and variable hierarchy
An empty GDS file can be created by createfn.gds():
library(gdsfmt)
gfile <- createfn.gds("test.gds")Now, a file handle associated with “test.gds” is saved in the R variable gfile.
The GDS file can contain a hierarchical structure to store multiple GDS variables (or GDS nodes) in the file, and various data types are allowed (see the document of add.gdsn()) including integer, floating-point number and character.
add.gdsn(gfile, "int", val=1:10000)
add.gdsn(gfile, "double", val=seq(1, 1000, 0.4))
add.gdsn(gfile, "character", val=c("int", "double", "logical", "factor"))
add.gdsn(gfile, "logical", val=rep(c(TRUE, FALSE, NA), 50), visible=FALSE)
add.gdsn(gfile, "factor", val=as.factor(c(NA, "AA", "CC")), visible=FALSE)
add.gdsn(gfile, "bit2", val=sample(0:3, 1000, replace=TRUE), storage="bit2")
# list and data.frame
add.gdsn(gfile, "list", val=list(X=1:10, Y=seq(1, 10, 0.25)))
add.gdsn(gfile, "data.frame", val=data.frame(X=1:19, Y=seq(1, 10, 0.5)))folder <- addfolder.gdsn(gfile, "folder")
add.gdsn(folder, "int", val=1:1000)
add.gdsn(folder, "double", val=seq(1, 100, 0.4), visible=FALSE)Users can display the file content by typing gfile or print(gfile):
gfile## File: /tmp/Rtmpzkev6e/Rbuild57c440e7b8d0/gdsfmt/vignettes/test.gds (1.2 KB)
## + [ ]
## |--+ int { Int32 10000, 40.0 KB }
## |--+ double { Float64 2498, 20.0 KB }
## |--+ character { VStr8 4, 26 bytes }
## |--+ bit2 { Bit2 1000, 250 bytes }
## |--+ list [ list ] *
## | |--+ X { Int32 10, 40 bytes }
## | |--+ Y { Float64 37, 296 bytes }
## |--+ data.frame [ data.frame ] *
## | |--+ X { Int32 19, 76 bytes }
## | |--+ Y { Float64 19, 152 bytes }
## |--+ folder [ ]
## | |--+ int { Int32 1000, 4.0 KB }print(gfile, ...) has an argument all to control the display of file content. By default, all=FALSE; if all=TRUE, to show all contents in the file including hidden variables or folders. The GDS variables logical, factor and folder/double are hidden.
print(gfile, all=TRUE)## File: /tmp/Rtmpzkev6e/Rbuild57c440e7b8d0/gdsfmt/vignettes/test.gds (68.9 KB)
## + [ ]
## |--+ int { Int32 10000, 40.0 KB }
## |--+ double { Float64 2498, 20.0 KB }
## |--+ character { VStr8 4, 26 bytes }
## |--+ logical { Int32,logical 150, 600 bytes } *
## |--+ factor { Int32,factor 3, 12 bytes } *
## |--+ bit2 { Bit2 1000, 250 bytes }
## |--+ list [ list ] *
## | |--+ X { Int32 10, 40 bytes }
## | |--+ Y { Float64 37, 296 bytes }
## |--+ data.frame [ data.frame ] *
## | |--+ X { Int32 19, 76 bytes }
## | |--+ Y { Float64 19, 152 bytes }
## |--+ folder [ ]
## | |--+ int { Int32 1000, 4.0 KB }
## | |--+ double { Float64 248, 2.0 KB } *The asterisk indicates attributes attached to a GDS variable. The attributes can be used in the R environment to interpret the variable as logical, factor, data.frame or list.
index.gdsn() can locate the GDS variable by a path:
index.gdsn(gfile, "int")## + int { Int32 10000, 40.0 KB }index.gdsn(gfile, "list/Y")## + list/Y { Float64 37, 296 bytes }index.gdsn(gfile, "folder/int")## + folder/int { Int32 1000, 4.0 KB }# close the GDS file
closefn.gds(gfile)Writing Data
Array-oriented data sets can be written to the GDS file. There are three possible ways to write data to a GDS variable.
gfile <- createfn.gds("test.gds")R function add.gdsn
Users could pass an R variable to the function add.gdsn() directly. read.gdsn() can read data from the GDS variable.
(n <- add.gdsn(gfile, "I1", val=matrix(1:15, nrow=3)))## + I1 { Int32 3x5, 60 bytes }read.gdsn(n)## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 4 7 10 13
## [2,] 2 5 8 11 14
## [3,] 3 6 9 12 15R function write.gdsn
Users can specify the arguments start and count to write a subset of data. -1 in count means the size of that dimension, and the corresponding element in start should be 1. The values in start and cound should be in the dimension range.
write.gdsn(n, rep(0,5), start=c(2,1), count=c(1,-1))
read.gdsn(n)## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 4 7 10 13
## [2,] 0 0 0 0 0
## [3,] 3 6 9 12 15R function append.gdsn
Users can append new data to an existing GDS variable.
append.gdsn(n, 16:24)
read.gdsn(n)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 1 4 7 10 13 16 19 22
## [2,] 0 0 0 0 0 17 20 23
## [3,] 3 6 9 12 15 18 21 24R function assign.gdsn
Users could call assign.gdsn() to replace specific values, subset or reorder the data variable.
# initialize
n <- add.gdsn(gfile, "mat", matrix(1:48, 6))
read.gdsn(n)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 1 7 13 19 25 31 37 43
## [2,] 2 8 14 20 26 32 38 44
## [3,] 3 9 15 21 27 33 39 45
## [4,] 4 10 16 22 28 34 40 46
## [5,] 5 11 17 23 29 35 41 47
## [6,] 6 12 18 24 30 36 42 48# substitute
assign.gdsn(n, .value=c(9:14,35:40), .substitute=NA)
read.gdsn(n)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 1 7 NA 19 25 31 NA 43
## [2,] 2 8 NA 20 26 32 NA 44
## [3,] 3 NA 15 21 27 33 NA 45
## [4,] 4 NA 16 22 28 34 NA 46
## [5,] 5 NA 17 23 29 NA 41 47
## [6,] 6 NA 18 24 30 NA 42 48# subset
assign.gdsn(n, seldim=list(rep(c(TRUE, FALSE),3), rep(c(FALSE, TRUE),4)))
read.gdsn(n)## [,1] [,2] [,3] [,4]
## [1,] 7 19 31 43
## [2,] NA 21 33 45
## [3,] NA 23 NA 47# initialize and subset
n <- add.gdsn(gfile, "mat", matrix(1:48, 6), replace=TRUE)
assign.gdsn(n, seldim=list(c(4,2,6,NA), c(5,6,NA,2,8,NA,4)))
read.gdsn(n)## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 28 34 NA 10 46 NA 22
## [2,] 26 32 NA 8 44 NA 20
## [3,] 30 36 NA 12 48 NA 24
## [4,] NA NA NA NA NA NA NA# initialize and sort into descending order
n <- add.gdsn(gfile, "mat", matrix(1:48, 6), replace=TRUE)
assign.gdsn(n, seldim=list(6:1, 8:1))
read.gdsn(n)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 48 42 36 30 24 18 12 6
## [2,] 47 41 35 29 23 17 11 5
## [3,] 46 40 34 28 22 16 10 4
## [4,] 45 39 33 27 21 15 9 3
## [5,] 44 38 32 26 20 14 8 2
## [6,] 43 37 31 25 19 13 7 1Create a large-scale data set
When the size of dataset is larger than the system memory, users can not add a GDS variable via add.gdsn() directly. If the dimension is pre-defined, users can specify the dimension size in add.gdsn() to allocate data space. Then call write.gdsn() to write a small subset of data space.
(n2 <- add.gdsn(gfile, "I2", storage="int", valdim=c(100, 2000)))## + I2 { Int32 100x2000, 800.0 KB }for (i in 1:2000)
{
write.gdsn(n2, seq.int(100*(i-1)+1, length.out=100),
start=c(1,i), count=c(-1,1))
}Or call append.gdsn() to append new data when the initial size is ZERO. If a compression algorithm is specified in add.gdsn() (e.g., compress=“ZIP”), users should call append.gdsn() instead of write.gdsn(), since data has to be compressed sequentially.
(n3 <- add.gdsn(gfile, "I3", storage="int", valdim=c(100, 0), compress="ZIP"))## + I3 { Int32 100x0 ZIP, 0 byte }for (i in 1:2000)
{
append.gdsn(n3, seq.int(100*(i-1)+1, length.out=100))
}
n3## + I3 { Int32 100x2000 ZIP(32.29%), 262.1 KB }# close the GDS file
closefn.gds(gfile)Reading Data
gfile <- createfn.gds("test.gds")
add.gdsn(gfile, "I1", val=matrix(1:20, nrow=4))
add.gdsn(gfile, "I2", val=1:100)
closefn.gds(gfile)read.gdsn() can load all data to an R variable in memory.
gfile <- openfn.gds("test.gds")
n <- index.gdsn(gfile, "I1")
read.gdsn(n)## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 5 9 13 17
## [2,] 2 6 10 14 18
## [3,] 3 7 11 15 19
## [4,] 4 8 12 16 20Subset reading read.gdsn and readex.gdsn
A subset of data can be specified via the arguments start and count in the R function read.gdsn. Or specify a list of logical vectors in readex.gdsn().
# read a subset
read.gdsn(n, start=c(2, 2), count=c(2, 3))## [,1] [,2] [,3]
## [1,] 6 10 14
## [2,] 7 11 15read.gdsn(n, start=c(2, 2), count=c(2, 3), .value=c(6,15), .substitute=NA)## [,1] [,2] [,3]
## [1,] NA 10 14
## [2,] 7 11 NA# read a subset
readex.gdsn(n, list(c(FALSE,TRUE,TRUE,FALSE), c(TRUE,FALSE,TRUE,FALSE,TRUE)))## [,1] [,2] [,3]
## [1,] 2 10 18
## [2,] 3 11 19readex.gdsn(n, list(c(1,4,3,NA), c(2,NA,3,1,3,1)))## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 5 NA 9 1 9 1
## [2,] 8 NA 12 4 12 4
## [3,] 7 NA 11 3 11 3
## [4,] NA NA NA NA NA NAreadex.gdsn(n, list(c(1,4,3,NA), c(2,NA,3,1,3,1)), .value=NA, .substitute=-1)## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 5 -1 9 1 9 1
## [2,] 8 -1 12 4 12 4
## [3,] 7 -1 11 3 11 3
## [4,] -1 -1 -1 -1 -1 -1Apply a user-defined function marginally
A user-defined function can be applied marginally to a GDS variable via apply.gdsn(). margin=1 indicates applying the function row by row, and margin=2 for applying the function column by column.
apply.gdsn(n, margin=1, FUN=print, as.is="none")## [1] 1 5 9 13 17
## [1] 2 6 10 14 18
## [1] 3 7 11 15 19
## [1] 4 8 12 16 20apply.gdsn(n, margin=2, FUN=print, as.is="none")## [1] 1 2 3 4
## [1] 5 6 7 8
## [1] 9 10 11 12
## [1] 13 14 15 16
## [1] 17 18 19 20# close the GDS file
closefn.gds(gfile)Examples
To create a simple GDS file,
gfile <- createfn.gds("test.gds")
n1 <- add.gdsn(gfile, "I1", val=1:100)
n2 <- add.gdsn(gfile, "I2", val=matrix(1:20, nrow=4))
gfile## File: /tmp/Rtmpzkev6e/Rbuild57c440e7b8d0/gdsfmt/vignettes/test.gds (224 bytes)
## + [ ]
## |--+ I1 { Int32 100, 400 bytes }
## |--+ I2 { Int32 4x5, 80 bytes }Output to a text file
apply.gdsn() can be used to export a GDS variable to a text file. If the GDS variable is a vector,
fout <- file("text.txt", "wt")
apply.gdsn(n1, 1, FUN=cat, as.is="none", file=fout, fill=TRUE)
close(fout)
scan("text.txt")## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
## [18] 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
## [35] 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51
## [52] 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68
## [69] 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85
## [86] 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100The arguments file and fill are defined in the function cat().
If the GDS variable is a matrix:
fout <- file("text.txt", "wt")
apply.gdsn(n2, 1, FUN=cat, as.is="none", file=fout, fill=4194304)
close(fout)
readLines("text.txt")## [1] "1 5 9 13 17" "2 6 10 14 18" "3 7 11 15 19" "4 8 12 16 20"The number 4194304 is the maximum number of columns on a line used in printing vectors.
Transpose a matrix
permdim.gdsn() can be used to transpose an array by permuting its dimensions. Or apply.gdsn() allows that the data returned from the user-defined function FUN is directly written to a target GDS node target.node, when as.is=“gdsnode” and target.node are both given. Little c in R is a generic function which combines its arguments, and it passes all data to the target GDS node in the following code:
n.t <- add.gdsn(gfile, "transpose", storage="int", valdim=c(5,0))
# apply the function over rows of matrix
apply.gdsn(n2, margin=1, FUN=c, as.is="gdsnode", target.node=n.t)
# matrix transpose
read.gdsn(n.t)## [,1] [,2] [,3] [,4]
## [1,] 1 2 3 4
## [2,] 5 6 7 8
## [3,] 9 10 11 12
## [4,] 13 14 15 16
## [5,] 17 18 19 20# close the GDS file
closefn.gds(gfile)Floating-point number vs. packed real number
In computing, floating point is a method of representing an approximation of a real number in a way that can support a trade-off between range and precision, which can be represented exactly is of the following form “significand \(\times\) 2exponent”. A packed real number in GDS format is defined as “int \(\times\) scale \(+\) offset”, where int can be a 8-bit, 16-bit or 32-bit signed interger. In some cases, the strategy of packed real numbers can significantly improve the compression ratio for real numbers.
set.seed(1000)
val <- sample(seq(0,1,0.001), 50000, replace=TRUE)
head(val)## [1] 0.328 0.759 0.114 0.691 0.516 0.067gfile <- createfn.gds("test.gds")
add.gdsn(gfile, "N1", val=val)
add.gdsn(gfile, "N2", val=val, compress="ZIP", closezip=TRUE)
add.gdsn(gfile, "N3", val=val, storage="float")
add.gdsn(gfile, "N4", val=val, storage="float", compress="ZIP", closezip=TRUE)
add.gdsn(gfile, "N5", val=val, storage="packedreal16", scale=0.001, offset=0)
add.gdsn(gfile, "N6", val=val, storage="packedreal16", scale=0.001, offset=0, compress="ZIP", closezip=TRUE)
gfile## File: /tmp/Rtmpzkev6e/Rbuild57c440e7b8d0/gdsfmt/vignettes/test.gds (793.7 KB)
## + [ ]
## |--+ N1 { Float64 50000, 400.0 KB }
## |--+ N2 { Float64 50000 ZIP(24.76%), 99.0 KB }
## |--+ N3 { Float32 50000, 200.0 KB }
## |--+ N4 { Float32 50000 ZIP(46.77%), 93.5 KB }
## |--+ N5 { PackedReal16 50000, 100.0 KB }
## |--+ N6 { PackedReal16 50000 ZIP(76.00%), 76.0 KB }| Variable | Type | Compression Method | Size | Ratio | Machine epsilon1 |
|---|---|---|---|---|---|
| N1 | 64-bit floating-point number | — | 400.0 KB | 100.0% | 0 |
| N2 | 64-bit floating-point number | ZIP | 99.0 KB | 24.8% | 0 |
| N3 | 32-bit floating-point number | — | 200.0 KB | 50.0% | 9.89e-09 |
| N4 | 32-bit floating-point number | ZIP | 93.5 KB | 23.4% | 9.89e-09 |
| N5 | 16-bit packed real number | — | 100.0 KB | 25.0% | 0 |
| N6 | 16-bit packed real number | ZIP | 76.0 KB | 19.0% | 0 |
1: the relative error due to rounding in floating point arithmetic.
# close the GDS file
closefn.gds(gfile)Limited random-access reading on compressed data
- 10,000,000 random 0,1 sequence of 32-bit integers
- in each 32 bits, one bit stores random 0,1 and others are ZERO
- lower bound of compression percentage is 1/32 = 3.125%
- Testing:
- of 10,000 random positions, read a 32-bit integer
- compression ratio is maximized for each method
- compression method: None, ZIP, ZIP_RA, LZ4 and LZ4_RA
- ZIP_RA and LZ4_RA: data are composed of multiple independent compressed blocks
set.seed(100)
# 10,000,000 random 0,1 sequence of 32-bit integers
val <- sample.int(2, 10*1000*1000, replace=TRUE) - 1L
table(val)## val
## 0 1
## 4999138 5000862# cteate a GDS file
f <- createfn.gds("test.gds")
# compression algorithms
compression <- c("", "ZIP.max", "ZIP_RA.max:16K", "LZ4.max", "LZ4_RA.max:16K")
# save
for (i in 1:length(compression))
print(add.gdsn(f, paste0("I", i), val=val, compress=compression[i], closezip=TRUE))
# close the file
closefn.gds(f)- System configuration:
- MacBook Pro, Retina, 13-inch, Late 2013, 2.8 GHz Intel Core i7, 16 GB 1600 MHz DDR3
# open the GDS file
f <- openfn.gds("test.gds")
# 10,000 random positions
set.seed(1000)
idx <- sample.int(length(val), 10000)
# enumerate each compression method
dat <- vector("list", length(compression))
for (i in seq_len(length(compression)))
{
cat("Compression:", compression[i], "\n")
n <- index.gdsn(f, paste0("I", i))
print(system.time({
dat[[i]] <- sapply(idx, FUN=function(k) read.gdsn(n, start=k, count=1L))
}))
}
# check
for (i in seq_len(length(compression)))
stopifnot(identical(dat[[i]], dat[[1L]]))
# close the file
closefn.gds(f)| Compression Method | Raw | ZIP | ZIP_RA | LZ4 | LZ4_RA |
|---|---|---|---|---|---|
| Data Size (MB) | 40.0 | 2.0 | 2.2 | 3.0 | 3.0 |
| Compression Percent | 100% | 5.08% | 5.42% | 7.39% | 7.60% |
| Reading Time (second) | 0.38 | 193.43 | 3.40 | 77.46 | 1.42 |
Stylish Terminal Output in R
If the R package crayon is installed in the R environment, print() can display the context of GDS file with different colours. For example, on Apple Mac, 
Users can disable crayon terminal output by options(gds.crayon=FALSE),
File: 1KG_autosome_phase3_shapeit2_mvncall_integrated_v5_20130502_genotypes.gds (3.4 GB)
+ [ ] *
|--+ sample.id { VStr8 2504 ZIP_RA(27.15%), 5.4 KB }
|--+ snp.id { Int32 81271745 ZIP_RA(34.58%), 112.4 MB }
|--+ snp.rs.id { VStr8 81271745 ZIP_RA(38.67%), 193.1 MB }
|--+ snp.position { Int32 81271745 ZIP_RA(39.73%), 129.1 MB }
|--+ snp.chromosome { VStr8 81271745 ZIP_RA(0.10%), 190.2 KB }
|--+ snp.allele { VStr8 81271745 ZIP_RA(17.05%), 57.3 MB }
|--+ genotype { Bit2 2504x81271745 ZIP_RA(5.66%), 2.9 GB } *
|--+ snp.annot [ ]
| |--+ qual { Float32 81271745 ZIP_RA(0.10%), 316.1 KB }
| |--+ filter { VStr8 81271745 ZIP_RA(0.15%), 592.0 KB }Session Information
sessionInfo()## R version 3.2.2 (2015-08-14)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 14.04.3 LTS
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] gdsfmt_1.6.2
##
## loaded via a namespace (and not attached):
## [1] magrittr_1.5 formatR_1.2.1 htmltools_0.2.6 tools_3.2.2
## [5] yaml_2.1.13 memoise_0.2.1 crayon_1.3.1 stringi_1.0-1
## [9] rmarkdown_0.8.1 knitr_1.11 stringr_1.0.0 digest_0.6.8
## [13] evaluate_0.8