# Define test file
#
suppressPackageStartupMessages({
    library(HDF5Array)
    library(data.table)
})
d <- '/s/project/patient_report/tidy_results/prokisch-full/savedObjects/Full_Prokisch_Analysis'
f <- 'zscore_psiSite.h5'
fp <- file.path(d,f)
h5ls(fp)
##   group           name       otype dclass           dim
## 0     / zscore_psiSite H5I_DATASET  FLOAT 3269955 x 146
da <- HDF5Array(fp, sub('\\.h5','',f))

Subsetting by rows

loading does not scale linearly

log10(nrow(da))
## [1] 6.514542
system.time(as.matrix(da[1:1e3,]))
##    user  system elapsed 
##   0.145   0.005   0.151
system.time(as.matrix(da[1:1e4,]))
##    user  system elapsed 
##   0.669   0.006   0.676
system.time(as.matrix(da[1:1e5,]))
##    user  system elapsed 
##   1.514   0.059   1.572

random access is slower than by blocks

system.time(as.matrix(da[sample(1:nrow(da), 1e3),]))
##    user  system elapsed 
##  14.531   0.284  14.815

Subsetting by columns

ncol(da)
## [1] 146
system.time(as.matrix(da[,1]))
##    user  system elapsed 
##   1.270   0.020   1.291
system.time(as.matrix(da[,1:10]))
##    user  system elapsed 
##   3.021   0.097   3.119
system.time(as.matrix(da[,1:25]))
##    user  system elapsed 
##   5.816   0.479   6.295
system.time(as.matrix(da[,1:50]))
##    user  system elapsed 
##   9.761   1.421  11.184

random access is slower than by blocks

system.time(as.matrix(da[,sample(1:ncol(da), 10)]))
##    user  system elapsed 
##   4.817   0.255   5.072
system.time(as.matrix(da[,sample(1:ncol(da), 25)]))
##    user  system elapsed 
##  11.094   0.539  11.632

Use case I: Particular rows for 1 column

system.time(as.matrix(da[,1]))
##    user  system elapsed 
##   1.283   0.037   1.320
system.time(as.matrix(da[sample(1:nrow(da), 1e3),1]))
##    user  system elapsed 
##   0.181   0.003   0.184
system.time(as.matrix(da[sample(1:nrow(da), 1e4),1]))
##    user  system elapsed 
##   7.999   0.002   8.009

If the row subset is faster, you could spend more time in a clever row selection. For us this could be no missing values, no Pvalues>0.5, …, to reduce the number of points we plot in the end. Typical size would be 5k to 10k points to plot.

Use case II: Particular rows for multiple columns

Extract multiple features (organized in different HDF5 files) for 1 patient (1 column per file), but only the relevant rows (random indices).

system.time(as.matrix(da[,1:10]))
##    user  system elapsed 
##    2.35    0.16    2.51
system.time(as.matrix(da[sample(1:nrow(da), 1e3), 1:10]))
##    user  system elapsed 
##   1.267   0.020   1.287
system.time(as.matrix(da[sample(1:nrow(da), 1e3), sample(1:ncol(da), 10)]))
##    user  system elapsed 
##   1.850   0.008   1.859
system.time(as.matrix(da[sample(1:nrow(da), 1e4), 1:10]))
##    user  system elapsed 
##  14.308   0.019  14.327
# very slow ~2min
# system.time(as.matrix(da[sample(1:nrow(da), 1e4), sample(1:ncol(da), 10)]))