| Title: | Read Genomic Data Files (VCF, BCF, BGEN, PGEN, BED, H5AD) into R/Rcpp in Chunks |
|---|---|
| Description: | Read genomic data files (i.e. VCF, BCF, BGEN, PGEN, BED, H5AD, DelayedArray) into R/Rcpp in chunks for analysis with Armadillo / Eigen / Rcpp libraries. Modern datasets are often too big to fit into memory, and many analyses operate on a small chunk features at a time. GenomicDataStream provides interfaces at both the C++ and R levels for out-of-core analysis. |
| Authors: | Gabriel Hoffman [aut, cre] (ORCID: <https://orcid.org/0000-0002-0957-0224>), Zilong Li [aut] (ORCID: <https://orcid.org/0000-0001-5859-2078>) |
| Maintainer: | Gabriel Hoffman <[email protected]> |
| License: | Artistic-2.0 |
| Version: | 0.99.28 |
| Built: | 2026-09-30 00:46:03 UTC |
| Source: | https://github.com/BiocStaging/GenomicDataStream |
Detected if end of stream is reached
atEndOfStream(x)atEndOfStream(x)
x |
|
if end of stream has been reached, return TRUE. Else FALSE
file <- system.file("extdata", "test.vcf.gz", package = "GenomicDataStream") # initialize obj <- GenomicDataStream(file, "DS", chunkSize = 5, initialize = TRUE) # loop until break while (1) { # get data chunk # data$X matrix with features as columns # data$info information about each feature as rows dat <- getNextChunk(obj) if (atEndOfStream(obj)) break print(dat$info) }file <- system.file("extdata", "test.vcf.gz", package = "GenomicDataStream") # initialize obj <- GenomicDataStream(file, "DS", chunkSize = 5, initialize = TRUE) # loop until break while (1) { # get data chunk # data$X matrix with features as columns # data$info information about each feature as rows dat <- getNextChunk(obj) if (atEndOfStream(obj)) break print(dat$info) }
Get genomic intervals chopping up regions into smaller chunks
chopChroms(x, nchunks = 10)chopChroms(x, nchunks = 10)
x |
|
nchunks |
number of chunks to create |
vector genome intervals
file <- system.file("extdata", "test.bed", package = "GenomicDataStream") # initialize obj <- GenomicDataStream(file, chunkSize = 5, initialize = TRUE) df <- getChromRanges(obj) chopChroms( df )file <- system.file("extdata", "test.bed", package = "GenomicDataStream") # initialize obj <- GenomicDataStream(file, chunkSize = 5, initialize = TRUE) df <- getChromRanges(obj) chopChroms( df )
Get array of chunk sizes
countChunks(x)countChunks(x)
x |
|
array of chunk sizes
file <- system.file("extdata", "test.vcf.gz", package = "GenomicDataStream") # initialize obj <- GenomicDataStream(file, "DS", chunkSize = 5, initialize = TRUE) countChunks(obj)file <- system.file("extdata", "test.vcf.gz", package = "GenomicDataStream") # initialize obj <- GenomicDataStream(file, "DS", chunkSize = 5, initialize = TRUE) countChunks(obj)
DashSVD for fast randomized SVD
dashSVD(x, k, p = 3, s = 20, rand = 1, byrow = FALSE, verbose = FALSE) ## S4 method for signature 'ANY' dashSVD(x, k, p = 3, s = 20, rand = 1, byrow = FALSE, verbose = FALSE)dashSVD(x, k, p = 3, s = 20, rand = 1, byrow = FALSE, verbose = FALSE) ## S4 method for signature 'ANY' dashSVD(x, k, p = 3, s = 20, rand = 1, byrow = FALSE, verbose = FALSE)
x |
data matrix |
k |
specifies the target rank of the low-rank decomposition. |
p |
number of additional power iterations (by default |
s |
oversampling parameter (by default |
rand |
randomization method. 1 indicates Gaussian sampling and 2 indicates uniform sampling |
byrow |
default |
verbose |
print messages |
PCA class storing d, u and v
Zilong Li [email protected], Gabriel Hoffman
Feng, X., Yu, W., Xie, Y. and Tang, J., 2024. Algorithm 1043: Faster randomized svd with dynamic shifts. ACM Transactions on Mathematical Software, 50(2), pp.1-27.
hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, `+`) } X <- hilbert(9)[, 1:6] dcmp1 <- svd(X) dcmp2 <- dashSVD(X, k=5) # compare singular values dcmp1$d dcmp2$dhilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, `+`) } X <- hilbert(9)[, 1:6] dcmp1 <- svd(X) dcmp2 <- dashSVD(X, k=5) # compare singular values dcmp1$d dcmp2$d
Get number of total features read from GenomicDataStream
featuresRead(x)featuresRead(x)
x |
|
total number of features read from the stream
file <- system.file("extdata", "test.vcf.gz", package = "GenomicDataStream") # initialize obj <- GenomicDataStream(file, "DS", chunkSize = 5, initialize = TRUE) # loop until break while (1) { # get data chunk # data$X matrix with features as columns # data$info information about each feature as rows dat <- getNextChunk(obj) if (atEndOfStream(obj)) break print(dat$info) } featuresRead(obj)file <- system.file("extdata", "test.vcf.gz", package = "GenomicDataStream") # initialize obj <- GenomicDataStream(file, "DS", chunkSize = 5, initialize = TRUE) # loop until break while (1) { # get data chunk # data$X matrix with features as columns # data$info information about each feature as rows dat <- getNextChunk(obj) if (atEndOfStream(obj)) break print(dat$info) } featuresRead(obj)
Read genomic data files (VCF, BCF, BGEN, h5ad) into R/Rcpp in chunks for analysis with Armadillo or Eigen libraries
Interface to GenomicDataStream C++ code
GenomicDataStream( file, field = "", region = "", samples = "-", MAF = 0, minVariance = 0, chunkSize = 10000, missingToMean = TRUE, initialize = FALSE )GenomicDataStream( file, field = "", region = "", samples = "-", MAF = 0, minVariance = 0, chunkSize = 10000, missingToMean = TRUE, initialize = FALSE )
file |
file in VCF/BCF/BGEN/PGEN format with index |
field |
field of VCF/BCF to read. Ignored for other file types |
region |
target in the format |
samples |
string of comma separated sample IDs to extract: |
MAF |
minor allele frequency filter applied to variants with max value <= 2. Use |
minVariance |
minimum variance filter applied to variants with max value > 2 |
chunkSize |
number of variants to return per chunk |
missingToMean |
if true, set missing values to the mean dosage value. if false, set to |
initialize |
default |
Variants are filtered using MAF if the max value is <=2, or minVariance otherwise
object of class GenomicDataStream
file <- system.file("extdata", "test.vcf.gz", package = "GenomicDataStream") # initialize obj <- GenomicDataStream(file, "DS", chunkSize = 5, initialize = TRUE) obj # loop until break while (1) { # get data chunk # data$X matrix with features as columns # data$info information about each feature as rows dat <- getNextChunk(obj) if (atEndOfStream(obj)) break print(dat$info) }file <- system.file("extdata", "test.vcf.gz", package = "GenomicDataStream") # initialize obj <- GenomicDataStream(file, "DS", chunkSize = 5, initialize = TRUE) obj # loop until break while (1) { # get data chunk # data$X matrix with features as columns # data$info information about each feature as rows dat <- getNextChunk(obj) if (atEndOfStream(obj)) break print(dat$info) }
Read genomic data files (VCF, BCF, BGEN, h5ad) into R/Rcpp in chunks for analysis with Armadillo or Eigen libraries
Get max and max postion for each chromosome
getChromRanges(x)getChromRanges(x)
x |
|
data.frame of chrom, start, end
file <- system.file("extdata", "test.bed", package = "GenomicDataStream") # initialize obj <- GenomicDataStream(file, chunkSize = 5, initialize = TRUE) getChromRanges( obj )file <- system.file("extdata", "test.bed", package = "GenomicDataStream") # initialize obj <- GenomicDataStream(file, chunkSize = 5, initialize = TRUE) getChromRanges( obj )
Get data chunk from GenomicDataStream
getNextChunk(x)getNextChunk(x)
x |
|
get data chunk as list with entries X storing a matrix with features as columns, and info storing information about each feature as rows
file <- system.file("extdata", "test.vcf.gz", package = "GenomicDataStream") # initialize obj <- GenomicDataStream(file, "DS", chunkSize = 5, initialize = TRUE) # loop until break while (1) { # get data chunk # data$X matrix with features as columns # data$info information about each feature as rows dat <- getNextChunk(obj) if (atEndOfStream(obj)) break print(dat$info) }file <- system.file("extdata", "test.vcf.gz", package = "GenomicDataStream") # initialize obj <- GenomicDataStream(file, "DS", chunkSize = 5, initialize = TRUE) # loop until break while (1) { # get data chunk # data$X matrix with features as columns # data$info information about each feature as rows dat <- getNextChunk(obj) if (atEndOfStream(obj)) break print(dat$info) }
BGEN uses sample order from the query, but VCF/BCF/PGEN uses order in file
getSampleNames(x)getSampleNames(x)
x |
|
Get sample names in order that the genotypes are extracted
array of string names
file <- system.file("extdata", "test.vcf.gz", package = "GenomicDataStream") obj <- GenomicDataStream(file, "DS", init=TRUE) getSampleNames(obj)file <- system.file("extdata", "test.vcf.gz", package = "GenomicDataStream") obj <- GenomicDataStream(file, "DS", init=TRUE) getSampleNames(obj)
Get location of each feature in GenomicDataStream
getVariantLocations(x)getVariantLocations(x)
x |
|
get data chunk as list with entries X storing a matrix with features as columns, and info storing information about each feature as rows
file <- system.file("extdata", "test.vcf.gz", package = "GenomicDataStream") # initialize obj <- GenomicDataStream(file, "DS", chunkSize = 5, initialize = TRUE) getVariantLocations(obj)file <- system.file("extdata", "test.vcf.gz", package = "GenomicDataStream") # initialize obj <- GenomicDataStream(file, "DS", chunkSize = 5, initialize = TRUE) getVariantLocations(obj)
Read file info from path to initialise stream. If already initialized, return to the beginning of the stream
initializeStream(x, region = NULL)initializeStream(x, region = NULL)
x |
|
region |
target in the format |
initialized GenomicDataStream
file <- system.file("extdata", "test.vcf.gz", package = "GenomicDataStream") obj <- GenomicDataStream(file, "DS", chunkSize = 5) # by default, GenomicDataStream is not initialized isInitialized(obj) # initialize obj <- initializeStream(obj) isInitialized(obj)file <- system.file("extdata", "test.vcf.gz", package = "GenomicDataStream") obj <- GenomicDataStream(file, "DS", chunkSize = 5) # by default, GenomicDataStream is not initialized isInitialized(obj) # initialize obj <- initializeStream(obj) isInitialized(obj)
If initialized, return TRUE, else FALSE
isInitialized(x)isInitialized(x)
x |
|
initialization status
file <- system.file("extdata", "test.vcf.gz", package = "GenomicDataStream") obj <- GenomicDataStream(file, "DS", chunkSize = 5) # by default, GenomicDataStream is not initialized isInitialized(obj) # initialize obj <- initializeStream(obj) isInitialized(obj)file <- system.file("extdata", "test.vcf.gz", package = "GenomicDataStream") obj <- GenomicDataStream(file, "DS", chunkSize = 5) # by default, GenomicDataStream is not initialized isInitialized(obj) # initialize obj <- initializeStream(obj) isInitialized(obj)
Modify signs of principal components so diagonal are always positive
normPC(U)normPC(U)
U |
matrix of principal components |
matrix of principal components with positive diagonal values
hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, `+`) } X <- hilbert(9)[, 1:6] k <- 4 dcmp <- svd(X, k, k) normPC( dcmp$u )hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, `+`) } X <- hilbert(9)[, 1:6] k <- 4 dcmp <- svd(X, k, k) normPC( dcmp$u )
Window-based Randomized SVD with PCAone
PCAone( x, k, ..., p = 7, s = 20, B = 64, threads = 4, threads2 = 4, scaleAndCenter = TRUE, shuffle = TRUE, verbose = TRUE ) ## S4 method for signature 'GenomicDataStream' PCAone( x, k, ..., p = 7, s = 20, B = 64, threads = 4, threads2 = 1, scaleAndCenter = TRUE, shuffle = TRUE, verbose = TRUE )PCAone( x, k, ..., p = 7, s = 20, B = 64, threads = 4, threads2 = 4, scaleAndCenter = TRUE, shuffle = TRUE, verbose = TRUE ) ## S4 method for signature 'GenomicDataStream' PCAone( x, k, ..., p = 7, s = 20, B = 64, threads = 4, threads2 = 1, scaleAndCenter = TRUE, shuffle = TRUE, verbose = TRUE )
x |
|
k |
integer; |
... |
other argument to control streaming |
p |
integer, optional; |
s |
integer, optional; |
B |
integer, optional; |
threads |
integer, optional; |
threads2 |
integer, optional; |
scaleAndCenter |
bool, optional; |
shuffle |
bool, optional; |
verbose |
string, optional; |
PCAone implements the window-based Randomized SVD proposed by Li, et al. (2023).
If scaleAndCenter, the data matrix is scaled so that each feature has a mean of zero and a cross-product of 1. In this case the sum of squares of all eigen values equals the number of features (i.e sum(d^2) = p).
Computational time is spent on two steps:
1) Reading data from disk and and processing data. Multiple chunks can be read and processed in parallel. This is controlled by setting threads
2) Updating PCA with current data chunk. Only one chunk can be processed at a time, but linear algebra operations can be parallelized. This is conrolled by setting threads2
If x is a matrix, PCA is performed using rows as features. Note: this is the _transpose_ of what svd() does.
PCAone returns a list containing the following three components:
array_like;
singular values; vector of length .
array_like;
left singular vectors; or dimensional array.
array_like;
right singular vectors; or dimensional array.
The singular vectors are not unique and only defined up to sign. If a left singular vector has its sign changed, changing the sign of the corresponding right vector gives an equivalent decomposition.
Zilong Li [email protected], Gabriel Hoffman
Li, Z., Meisner, J., & Albrechtsen, A. (2023). Fast and accurate out-of-core PCA framework for large scale biobank data. Genome Research, 33(9), 1599-1608. doi:10.1101/gr.277525.122.
# Example: PCA on genotype data file <- system.file("extdata", "test.vcf.gz", package = "GenomicDataStream") obj <- GenomicDataStream(file, "DS", chunkSize = 3) res <- PCAone(obj, k=5, threads=1) res par(pty="s") plot(res)# Example: PCA on genotype data file <- system.file("extdata", "test.vcf.gz", package = "GenomicDataStream") obj <- GenomicDataStream(file, "DS", chunkSize = 3) res <- PCAone(obj, k=5, threads=1) res par(pty="s") plot(res)
Evaluate performance of PC estimates compared to true PC values
perfMetric(U, U_est, k = ncol(U_est), metric = c("MEV", "minSSE"))perfMetric(U, U_est, k = ncol(U_est), metric = c("MEV", "minSSE"))
U |
true eigen values |
U_est |
estimated eigen-values |
k |
number of PCs to evaluate |
metric |
evaluate the accuracy of the estimated PCs compared to the true PCs using mean explained variance ( |
See performance metrics described by Li, et al. (2023)
performance metric
Li, Z., Meisner, J., & Albrechtsen, A. (2023). Fast and accurate out-of-core PCA framework for large scale biobank data. Genome Research, 33(9), 1599-1608. doi:10.1101/gr.277525.122.
hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, `+`) } X <- hilbert(9)[, 1:6] k <- 4 # Compute SVD using two methods dcmp1 <- svd(X, k, k) dcmp2 <- dashSVD( X, k=k) # Mean variance explained is 1 perfMetric(dcmp1$u, dcmp2$u, metric = "MEV") # minimum of sum of squared errors zero perfMetric(dcmp1$u, dcmp2$u, metric = "minSSE")hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, `+`) } X <- hilbert(9)[, 1:6] k <- 4 # Compute SVD using two methods dcmp1 <- svd(X, k, k) dcmp2 <- dashSVD( X, k=k) # Mean variance explained is 1 perfMetric(dcmp1$u, dcmp2$u, metric = "MEV") # minimum of sum of squared errors zero perfMetric(dcmp1$u, dcmp2$u, metric = "minSSE")
Plot PCAone
## S4 method for signature 'PCA,ANY' plot(x, ..., main = "scree plot")## S4 method for signature 'PCA,ANY' plot(x, ..., main = "scree plot")
x |
|
... |
other arguments |
main |
title |
plot
# Example: PCA on genotype data file <- system.file("extdata", "test.vcf.gz", package = "GenomicDataStream") obj <- GenomicDataStream(file, "DS", chunkSize = 3) res <- PCAone(obj, k=5, threads=1) plot(res)# Example: PCA on genotype data file <- system.file("extdata", "test.vcf.gz", package = "GenomicDataStream") obj <- GenomicDataStream(file, "DS", chunkSize = 3) res <- PCAone(obj, k=5, threads=1) plot(res)
Print object
Print object
## S4 method for signature 'GenomicDataStream' print(x, ...) ## S4 method for signature 'PCA' print(x, ...)## S4 method for signature 'GenomicDataStream' print(x, ...) ## S4 method for signature 'PCA' print(x, ...)
x |
|
... |
other arguments |
printed strings
Read H5AD as SingleCellExperiment where counts is a file-backed DelayedArray
readH5AD(file, layer = NULL, ondisk = TRUE, verbose = FALSE, raw = FALSE)readH5AD(file, layer = NULL, ondisk = TRUE, verbose = FALSE, raw = FALSE)
file |
H5AD file |
layer |
|
ondisk |
if |
verbose |
print messages |
raw |
if |
Uses HDF5Array::H5ADMatrix() to read counts as a file-backed DelayedArray, and anndataR::read_h5ad() to read all other data from H5AD.
SingleCellExperiment
file <- system.file("extdata", "example.h5ad", package = "anndataR") sce <- readH5AD(file)file <- system.file("extdata", "example.h5ad", package = "anndataR") sce <- readH5AD(file)
BGEN uses sample order from the query, but VCF/BCF/PGEN uses order in file
## S4 method for signature 'GenomicDataStream' rownames(x, do.NULL = TRUE, prefix = "row")## S4 method for signature 'GenomicDataStream' rownames(x, do.NULL = TRUE, prefix = "row")
x |
|
do.NULL |
not used |
prefix |
not used |
Get rownames (i.e. sample names) in order that the samples are extracted
array of string names
file <- system.file("extdata", "test.vcf.gz", package = "GenomicDataStream") obj <- GenomicDataStream(file, "DS", init=TRUE) rownames(obj)file <- system.file("extdata", "test.vcf.gz", package = "GenomicDataStream") obj <- GenomicDataStream(file, "DS", init=TRUE) rownames(obj)
Set chunk size for existing GenomicDataStream
setChunkSize(x, chunkSize)setChunkSize(x, chunkSize)
x |
|
chunkSize |
positive integer |
none
file <- system.file("extdata", "test.vcf.gz", package = "GenomicDataStream") obj <- GenomicDataStream(file, "DS", chunkSize = 5, initialize=TRUE) obj <- setChunkSize(obj, 200)file <- system.file("extdata", "test.vcf.gz", package = "GenomicDataStream") obj <- GenomicDataStream(file, "DS", chunkSize = 5, initialize=TRUE) obj <- setChunkSize(obj, 200)
If GenomicDataStream is already initialized, set new query in C++ backend. Otherwise substitute region values
setRegion(x, region)setRegion(x, region)
x |
|
region |
target in the format |
Set regions of GenomicDataStream
GenomicDataStream with region set
file <- system.file("extdata", "test.vcf.gz", package = "GenomicDataStream") obj <- GenomicDataStream(file, "DS", chunkSize = 5) # by default, GenomicDataStream is not initialized setRegion(obj, "1:10000-12000")file <- system.file("extdata", "test.vcf.gz", package = "GenomicDataStream") obj <- GenomicDataStream(file, "DS", chunkSize = 5) # by default, GenomicDataStream is not initialized setRegion(obj, "1:10000-12000")
Show object
Show object
## S4 method for signature 'GenomicDataStream' show(object) ## S4 method for signature 'PCA' show(object)## S4 method for signature 'GenomicDataStream' show(object) ## S4 method for signature 'PCA' show(object)
object |
|
printed strings
Standardize mean and variance of matrix columns in place
standardize_in_place(X, center = TRUE, scale = TRUE)standardize_in_place(X, center = TRUE, scale = TRUE)
X |
matrix |
center |
boolean, TRUE indices center columns |
scale |
boolean, TRUE indices scale columns |
none, matrix is standardized in place
X <- matrix(runif(12), 3, 4) standardize_in_place(X) scale(X)X <- matrix(runif(12), 3, 4) standardize_in_place(X) scale(X)
Read through stream to get size of each chunk, IDs of variants and samples
summarizeChunks(x, threads = 4)summarizeChunks(x, threads = 4)
x |
|
threads |
number of threads |
list of info about stream
file <- system.file("extdata", "test.vcf.gz", package = "GenomicDataStream") # initialize obj <- GenomicDataStream(file, "DS", chunkSize = 3, initialize = TRUE) summarizeChunks(obj)file <- system.file("extdata", "test.vcf.gz", package = "GenomicDataStream") # initialize obj <- GenomicDataStream(file, "DS", chunkSize = 3, initialize = TRUE) summarizeChunks(obj)