First of all we need to install NewWave:
NewWave is a new package that assumes a Negative Binomial distributions for dimensionality reduction and batch effect removal. In order to reduce the memory consumption it uses a PSOCK cluster combined with the R package SharedObject that allow to share a matrix between different cores without memory duplication. Thanks to that we can massively parallelize the estimation process with huge benefit in terms of time consumption. We can reduce even more the time consumption using some minibatch approaches on the different steps of the optimization.
I am going to show how to use NewWave with example data generated with Splatter.
params <- newSplatParams()
N=500
set.seed(1234)
data <- splatSimulateGroups(params,batchCells=c(N/2,N/2),
group.prob = rep(0.1,10),
de.prob = 0.2,
verbose = FALSE)
Now we have a dataset with 500 cells and 10000 genes, I will use only the 500 most variable genes. NewWave takes as input raw data, not normalized.
set.seed(12359)
hvg <- rowVars(counts(data))
names(hvg) <- rownames(counts(data))
data <- data[names(sort(hvg,decreasing=TRUE))[1:500],]
As you can see there is a variable called batch in the colData section.
colData(data)
#> DataFrame with 500 rows and 4 columns
#> Cell Batch Group ExpLibSize
#> <character> <character> <factor> <numeric>
#> Cell1 Cell1 Batch1 Group8 63779.0
#> Cell2 Cell2 Batch1 Group5 64831.8
#> Cell3 Cell3 Batch1 Group4 62138.0
#> Cell4 Cell4 Batch1 Group2 49935.3
#> Cell5 Cell5 Batch1 Group1 68815.5
#> ... ... ... ... ...
#> Cell496 Cell496 Batch2 Group7 62672.0
#> Cell497 Cell497 Batch2 Group1 67586.5
#> Cell498 Cell498 Batch2 Group1 65997.7
#> Cell499 Cell499 Batch2 Group1 63210.5
#> Cell500 Cell500 Batch2 Group10 78498.3
IMPORTANT: For batch effecr removal the batch variable must be a factor
We also have a variable called Group that represent the cell type labels.
We can see the how the cells are distributed between group and batch
There is a clear batch effect between the cells.
Let’s try to correct it.
I am going to show different implementation and the suggested way to use them with the given hardware.
Some advise:
This is the way to insert the batch variable, in the same manner can be inserted other cell-related variable and if you need some gene related variable those can be inserted in V.
res <- newWave(data,X = "~Batch", K=10, verbose = TRUE)
#> Time of setup
#> user system elapsed
#> 0.021 0.000 0.701
#> Time of initialization
#> user system elapsed
#> 0.143 0.031 1.233
#> Iteration 1
#> penalized log-likelihood = -1289835.22540995
#> Time of dispersion optimization
#> user system elapsed
#> 0.818 0.021 0.839
#> after optimize dispersion = -1054268.10851508
#> Time of right optimization
#> user system elapsed
#> 0.002 0.000 12.037
#> after right optimization= -1053518.14207601
#> after orthogonalization = -1053518.12246619
#> Time of left optimization
#> user system elapsed
#> 0.002 0.001 11.889
#> after left optimization= -1053223.79948399
#> after orthogonalization = -1053223.79705115
#> Iteration 2
#> penalized log-likelihood = -1053223.79705115
#> Time of dispersion optimization
#> user system elapsed
#> 1.167 0.020 1.188
#> after optimize dispersion = -1053217.41886092
#> Time of right optimization
#> user system elapsed
#> 0.002 0.000 7.711
#> after right optimization= -1053182.80920191
#> after orthogonalization = -1053182.8056879
#> Time of left optimization
#> user system elapsed
#> 0.001 0.000 5.944
#> after left optimization= -1053168.85874376
#> after orthogonalization = -1053168.85866781
In order to make it faster you can increase the number of cores using “children” parameter:
res2 <- newWave(data,X = "~Batch", K=10, verbose = TRUE, children=2)
#> Time of setup
#> user system elapsed
#> 0.015 0.005 0.656
#> Time of initialization
#> user system elapsed
#> 0.069 0.008 0.685
#> Iteration 1
#> penalized log-likelihood = -1289835.2252412
#> Time of dispersion optimization
#> user system elapsed
#> 1.086 0.035 1.122
#> after optimize dispersion = -1054268.10715249
#> Time of right optimization
#> user system elapsed
#> 0.000 0.002 6.612
#> after right optimization= -1053518.14580683
#> after orthogonalization = -1053518.12621059
#> Time of left optimization
#> user system elapsed
#> 0.001 0.001 3.431
#> after left optimization= -1053223.81480013
#> after orthogonalization = -1053223.81237923
#> Iteration 2
#> penalized log-likelihood = -1053223.81237923
#> Time of dispersion optimization
#> user system elapsed
#> 0.676 0.004 0.680
#> after optimize dispersion = -1053217.43421444
#> Time of right optimization
#> user system elapsed
#> 0.002 0.000 3.045
#> after right optimization= -1053182.77273083
#> after orthogonalization = -1053182.76920975
#> Time of left optimization
#> user system elapsed
#> 0.002 0.000 3.982
#> after left optimization= -1053168.8829841
#> after orthogonalization = -1053168.88291862
If you do not have an high number of cores to run newWave this is the fastest way to run. The optimization process is done by three process itereated until convercence.
Each of these three steps can be accelerated using mini batch, the number of observation is settled with these parameters:
res3 <- newWave(data,X = "~Batch", verbose = TRUE,K=10, children=2,
n_gene_disp = 100, n_gene_par = 100, n_cell_par = 100)
#> Time of setup
#> user system elapsed
#> 0.017 0.000 0.571
#> Time of initialization
#> user system elapsed
#> 0.070 0.000 0.589
#> Iteration 1
#> penalized log-likelihood = -1289835.22547275
#> Time of dispersion optimization
#> user system elapsed
#> 1.106 0.008 1.115
#> after optimize dispersion = -1054268.11127105
#> Time of right optimization
#> user system elapsed
#> 0.003 0.000 6.667
#> after right optimization= -1053518.1672765
#> after orthogonalization = -1053518.14767281
#> Time of left optimization
#> user system elapsed
#> 0.002 0.000 6.000
#> after left optimization= -1053223.82570555
#> after orthogonalization = -1053223.82326886
#> Iteration 2
#> penalized log-likelihood = -1053223.82326886
#> Time of dispersion optimization
#> user system elapsed
#> 0.254 0.008 0.262
#> after optimize dispersion = -1053223.82326886
#> Time of right optimization
#> user system elapsed
#> 0.001 0.000 0.679
#> after right optimization= -1053218.54190017
#> after orthogonalization = -1053218.54151242
#> Time of left optimization
#> user system elapsed
#> 0.002 0.000 0.349
#> after left optimization= -1053218.34201849
#> after orthogonalization = -1053218.34199003
If you have a lot of core disposable or you want to estimate a genewise dispersion parameter this is the fastes configuration:
res3 <- newWave(data,X = "~Batch", verbose = TRUE,K=10, children=2,
n_gene_par = 100, n_cell_par = 100, commondispersion = FALSE)
#> Time of setup
#> user system elapsed
#> 0.019 0.000 0.383
#> Time of initialization
#> user system elapsed
#> 0.080 0.004 0.462
#> Iteration 1
#> penalized log-likelihood = -1289835.2254766
#> Time of dispersion optimization
#> user system elapsed
#> 0.635 0.032 0.667
#> after optimize dispersion = -1054268.11167594
#> Time of right optimization
#> user system elapsed
#> 0.002 0.001 3.876
#> after right optimization= -1053518.14648746
#> after orthogonalization = -1053518.12688126
#> Time of left optimization
#> user system elapsed
#> 0.001 0.000 3.271
#> after left optimization= -1053223.79955525
#> after orthogonalization = -1053223.79713501
#> Iteration 2
#> penalized log-likelihood = -1053223.79713501
#> Time of dispersion optimization
#> user system elapsed
#> 0.072 0.008 0.538
#> after optimize dispersion = -1049102.25556613
#> Time of right optimization
#> user system elapsed
#> 0.001 0.000 0.618
#> after right optimization= -1049091.00200616
#> after orthogonalization = -1049090.9854972
#> Time of left optimization
#> user system elapsed
#> 0.002 0.000 0.679
#> after left optimization= -1049055.57694986
#> after orthogonalization = -1049055.57595901
#> Iteration 3
#> penalized log-likelihood = -1049055.57595901
#> Time of dispersion optimization
#> user system elapsed
#> 0.047 0.008 0.220
#> after optimize dispersion = -1049055.58898934
#> Time of right optimization
#> user system elapsed
#> 0.001 0.000 0.595
#> after right optimization= -1049050.42200897
#> after orthogonalization = -1049050.42164628
#> Time of left optimization
#> user system elapsed
#> 0.001 0.001 0.676
#> after left optimization= -1049018.65993138
#> after orthogonalization = -1049018.65906618
NB: do not use n_gene_disp in this case, it will slower the computation.
Now I can use the latent dimension rapresentation for visualization purpose:
latent <- reducedDim(res)
tsne_latent <- data.frame(Rtsne(latent)$Y)
tsne_latent$batch <- data$Batch
tsne_latent$group <- data$Group
or for clustering:
sessionInfo()
#> R version 4.4.0 RC (2024-04-16 r86468)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 22.04.4 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.20-bioc/R/lib/libRblas.so
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_GB 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
#>
#> time zone: America/New_York
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] NewWave_1.15.0 mclust_6.1.1
#> [3] ggplot2_3.5.1 Rtsne_0.17
#> [5] irlba_2.3.5.1 Matrix_1.7-0
#> [7] splatter_1.29.0 SingleCellExperiment_1.27.0
#> [9] SummarizedExperiment_1.35.0 Biobase_2.65.0
#> [11] GenomicRanges_1.57.0 GenomeInfoDb_1.41.0
#> [13] IRanges_2.39.0 S4Vectors_0.43.0
#> [15] BiocGenerics_0.51.0 MatrixGenerics_1.17.0
#> [17] matrixStats_1.3.0
#>
#> loaded via a namespace (and not attached):
#> [1] gtable_0.3.5 xfun_0.43 bslib_0.7.0
#> [4] lattice_0.22-6 vctrs_0.6.5 tools_4.4.0
#> [7] generics_0.1.3 parallel_4.4.0 tibble_3.2.1
#> [10] fansi_1.0.6 highr_0.10 pkgconfig_2.0.3
#> [13] SharedObject_1.19.0 checkmate_2.3.1 lifecycle_1.0.4
#> [16] GenomeInfoDbData_1.2.12 farver_2.1.1 compiler_4.4.0
#> [19] munsell_0.5.1 codetools_0.2-20 htmltools_0.5.8.1
#> [22] sass_0.4.9 yaml_2.3.8 pillar_1.9.0
#> [25] crayon_1.5.2 jquerylib_0.1.4 BiocParallel_1.39.0
#> [28] DelayedArray_0.31.0 cachem_1.0.8 abind_1.4-5
#> [31] rsvd_1.0.5 tidyselect_1.2.1 locfit_1.5-9.9
#> [34] digest_0.6.35 BiocSingular_1.21.0 dplyr_1.1.4
#> [37] labeling_0.4.3 fastmap_1.1.1 grid_4.4.0
#> [40] colorspace_2.1-0 cli_3.6.2 SparseArray_1.5.0
#> [43] magrittr_2.0.3 S4Arrays_1.5.0 utf8_1.2.4
#> [46] withr_3.0.0 UCSC.utils_1.1.0 scales_1.3.0
#> [49] backports_1.4.1 rmarkdown_2.26 XVector_0.45.0
#> [52] httr_1.4.7 beachmat_2.21.0 ScaledMatrix_1.13.0
#> [55] evaluate_0.23 knitr_1.46 rlang_1.1.3
#> [58] Rcpp_1.0.12 glue_1.7.0 jsonlite_1.8.8
#> [61] R6_2.5.1 zlibbioc_1.51.0