# Large-scale single-cell RNA-seq data manipulation with GDS files

#### Dr. Xiuwen Zheng (Genomics Research Center, AbbVie)

#### Dec 2020

## Introduction

The SCArray package provides large-scale single-cell RNA-seq data manipulation using Genomic Data Structure (GDS) files. It combines dense/sparse matrices stored in GDS files and the Bioconductor infrastructure framework (SingleCellExperiment and DelayedArray) to provide out-of-memory data storage and manipulation using the R programming language. As shown in the figure, SCArray provides a `SingleCellExperiment` object for downstream data analyses. GDS is an alternative to HDF5. Unlike HDF5, GDS supports the direct storage of a sparse matrix without converting it to multiple vectors.

Figure 1: Workflow of SCArray

## Installation

- Requires R (>= v3.5.0), gdsfmt (>= v1.24.0)

- Bioconductor repository

To install this package, start R and enter:

    if (!requireNamespace("BiocManager", quietly=TRUE))
        install.packages("BiocManager")
    BiocManager::install("SCArray")

## Format conversion

### Conversion from SingleCellExperiment

The SCArray package can convert a single-cell experiment object (SingleCellExperiment) to a GDS file using the function `scConvGDS()`. For example,

    suppressPackageStartupMessages(library(SCArray))
    suppressPackageStartupMessages(library(SingleCellExperiment))

    # load a SingleCellExperiment object
    fn <- system.file("extdata", "LaMannoBrainSub.rds", package="SCArray")
    sce <- readRDS(fn)

    # convert to a GDS file
    scConvGDS(sce, "test.gds")

    ## Output: test.gds
    ## Compression: LZMA_RA
    ## Dimension: 18219 x 25
    ## Assay List:
    ##     counts  |+ counts   { SparseReal32 18219x25 LZMA_ra(13.1%), 75.6K }
    ## rowData:
    ## colData:
    ##     CELL_ID
    ##     Cell_type
    ## Done.
    ## Clean up the fragments of GDS file:
    ##     open the file 'test.gds' (138.1K)
    ##     # of fragments: 29
    ##     save to 'test.gds.tmp'
    ##     rename 'test.gds.tmp' (137.9K, reduced: 168B)
    ##     # of fragments: 15

    # list data structure in the GDS file
    (f <- scOpen("test.gds")); scClose(f)

    ## Object of class "SCArrayFileClass"
    ## File: C:\Users\biocbuild\bbs-3.14-bioc\tmpdir\RtmpqAfIGG\Rbuild1f60125d6ea9\SCArray\vignettes\test.gds (137.9K)
    ## +    [  ] *
    ## |--+ feature.id   { Str8 18219 LZMA_ra(48.6%), 60.5K }
    ## |--+ sample.id   { Str8 25 LZMA_ra(40.0%), 157B }
    ## |--+ counts   { SparseReal32 18219x25 LZMA_ra(13.1%), 75.6K }
    ## |--+ feature.data   [  ]
    ## |--+ sample.data   [  ]
    ## |  |--+ CELL_ID   { Str8 25 LZMA_ra(40.0%), 157B }
    ## |  \--+ Cell_type   { Str8 25 LZMA_ra(64.3%), 133B }
    ## \--+ meta.data   [  ]

### Conversion from a matrix

The input of `scConvGDS()` can be a dense or sparse matrix for count data:

    library(Matrix)

    cnt <- matrix(0, nrow=4, ncol=8)
    set.seed(100); cnt[sample.int(length(cnt), 8)] <- rpois(8, 4)
    (cnt <- as(cnt, "dgCMatrix"))

    ## 4 x 8 sparse Matrix of class "dgCMatrix"
    ##
    ## [1,] . . . . . . . 6
    ## [2,] 3 1 . . . 4 . .
    ## [3,] . . . . . 3 . 4
    ## [4,] 4 . 3 . . . . .

    # convert to a GDS file
    scConvGDS(cnt, "test.gds")

    ## Output: test.gds
    ## Compression: LZMA_RA
    ## Dimension: 4 x 8
    ## Assay List:
    ##     counts  |+ counts   { SparseReal32 4x8 LZMA_ra(159.4%), 109B }
    ## Done.
    ## Clean up the fragments of GDS file:
    ##     open the file 'test.gds' (1.4K)
    ##     # of fragments: 17
    ##     save to 'test.gds.tmp'
    ##     rename 'test.gds.tmp' (1.3K, reduced: 84B)
    ##     # of fragments: 10

## Single cell datasets

When a single-cell GDS file is available, users can use `scExperiment()` to load a SingleCellExperiment object from the GDS file. The assay data in the SingleCellExperiment object are DelayedMatrix objects to avoid the memory limit.

    # a GDS file in the SCArray package
    (fn <- system.file("extdata", "LaMannoBrainData.gds", package="SCArray"))

    ## [1] "C:/Users/biocbuild/bbs-3.14-bioc/tmpdir/RtmpqAfIGG/Rinst1f6076a04956/SCArray/extdata/LaMannoBrainData.gds"

    # load a SingleCellExperiment object from the file
    sce <- scExperiment(fn)
    sce

    ## class: SingleCellExperiment
    ## dim: 12000 243
    ## metadata(0):
    ## assays(1): counts
    ## rownames(12000): Rp1 Sox17 ... Efhd2 Fhad1
    ## rowData names(0):
    ## colnames(243): 1772072122_A04 1772072122_A05 ... 1772099011_H05
    ##   1772099012_E04
    ## colData names(2): CELL_ID Cell_type
    ## reducedDimNames(0):
    ## mainExpName: NULL
    ## altExpNames(0):

    # it is a DelayedMatrix (the whole matrix is not loaded)
    assays(sce)$counts

    ## <12000 x 243> sparse matrix of class SC_GDSMatrix and type "double":
    ##         1772072122_A04 1772072122_A05 ... 1772099011_H05 1772099012_E04
    ##     Rp1              0              0   .              0              0
    ##   Sox17              0              0   .              0              0
    ##  Mrpl15              1              2   .              2              2
    ##  Lypla1              0              0   .              0              1
    ##   Tcea1              1              0   .              6              1
    ##     ...              .              .   .              .              .
    ##   Agmat              0              0   .              0              0
    ## Dnajc16              0              0   .              0              0
    ##   Casp9              0              0   .              0              0
    ##   Efhd2              0              0   .              1              1
    ##   Fhad1              0              0   .              1              0

    # column data
    colData(sce)

    ## DataFrame with 243 rows and 2 columns
    ##                       CELL_ID   Cell_type
    ##                   <character> <character>
    ## 1772072122_A04 1772072122_A04     DA-VTA4
    ## 1772072122_A05 1772072122_A05     DA-VTA2
    ## 1772072122_A06 1772072122_A06     DA-VTA3
    ## 1772072122_A07 1772072122_A07      DA-SNC
    ## 1772072122_A08 1772072122_A08     DA-VTA4
    ## ...                       ...         ...
    ## 1772099011_D01 1772099011_D01     DA-VTA1
    ## 1772099011_F04 1772099011_F04     DA-VTA2
    ## 1772099011_G07 1772099011_G07     DA-VTA4
    ## 1772099011_H05 1772099011_H05      DA-SNC
    ## 1772099012_E04 1772099012_E04     DA-VTA4

    # row data
    rowData(sce)

    ## DataFrame with 12000 rows and 0 columns

## Data Manipulation and Analysis

SCArray provides a `SingleCellExperiment` object for downstream data analyses. At first, we create a log count matrix `logcnt` from the count matrix. Note that `logcnt` is also a DelayedMatrix without actually generating the whole matrix.

    cnt <- assays(sce)$counts
    logcnt <- log2(cnt + 1)
    assays(sce)$logcounts <- logcnt
    logcnt

    ## <12000 x 243> sparse matrix of class DelayedMatrix and type "double":
    ##         1772072122_A04 1772072122_A05 ... 1772099011_H05 1772099012_E04
    ##     Rp1       0.000000       0.000000   .       0.000000       0.000000
    ##   Sox17       0.000000       0.000000   .       0.000000       0.000000
    ##  Mrpl15       1.000000       1.584963   .       1.584963       1.584963
    ##  Lypla1       0.000000       0.000000   .       0.000000       1.000000
    ##   Tcea1       1.000000       0.000000   .       2.807355       1.000000
    ##     ...              .              .   .              .              .
    ##   Agmat              0              0   .              0              0
    ## Dnajc16              0              0   .              0              0
    ##   Casp9              0              0   .              0              0
    ##   Efhd2              0              0   .              1              1
    ##   Fhad1              0              0   .              1              0

### Mean for each column or row

The DelayedMatrixStats package provides functions operating on rows and columns of DelayedMatrix objects. For example, we can calculate the mean for each column or row of the log count matrix.

    suppressPackageStartupMessages(library(DelayedMatrixStats))

    col_mean <- DelayedMatrixStats::colMeans2(logcnt)
    str(col_mean)

    ##  num [1:243] 0.261 0.138 0.238 0.259 0.143 ...

    row_mean <- DelayedMatrixStats::rowMeans2(logcnt)
    str(row_mean)

    ##  num [1:12000] 0 0.00652 0.81827 0.47055 1.33912 ...

### UMAP analysis

The scater package can perform the uniform manifold approximation and projection (UMAP) for the cell data, based on the data in a SingleCellExperiment object.

    suppressPackageStartupMessages(library(scater))

    # run umap analysis
    sce <- runUMAP(sce)

`plotReducedDim()` plots cell-level reduced dimension results (UMAP) stored in the SingleCellExperiment object:

    plotReducedDim(sce, dimred="UMAP")

## Session Info

    # print version information about R, the OS and attached or loaded packages
    sessionInfo()

    ## R version 4.1.2 (2021-11-01)
    ## Platform: x86_64-w64-mingw32/x64 (64-bit)
    ## Running under: Windows Server 2012 R2 x64 (build 9600)
    ##
    ## Matrix products: default
    ##
    ## locale:
    ## [1] LC_COLLATE=C
    ## [2] LC_CTYPE=English_United States.1252
    ## [3] LC_MONETARY=English_United States.1252
    ## [4] LC_NUMERIC=C
    ## [5] LC_TIME=English_United States.1252
    ##
    ## attached base packages:
    ## [1] stats4    stats     graphics  grDevices utils     datasets  methods
    ## [8] base
    ##
    ## other attached packages:
    ##  [1] scater_1.22.0               ggplot2_3.3.5
    ##  [3] scuttle_1.4.0               DelayedMatrixStats_1.16.0
    ##  [5] SingleCellExperiment_1.16.0 SummarizedExperiment_1.24.0
    ##  [7] Biobase_2.54.0              GenomicRanges_1.46.1
    ##  [9] GenomeInfoDb_1.30.0         SCArray_1.2.1
    ## [11] DelayedArray_0.20.0         IRanges_2.28.0
    ## [13] S4Vectors_0.32.3            MatrixGenerics_1.6.0
    ## [15] matrixStats_0.61.0          BiocGenerics_0.40.0
    ## [17] Matrix_1.4-0                gdsfmt_1.30.0
    ##
    ## loaded via a namespace (and not attached):
    ##  [1] viridis_0.6.2           sass_0.4.0              BiocSingular_1.10.0
    ##  [4] viridisLite_0.4.0       jsonlite_1.7.3          bslib_0.3.1
    ##  [7] assertthat_0.2.1        highr_0.9               GenomeInfoDbData_1.2.7
    ## [10] vipor_0.4.5             yaml_2.2.2              ggrepel_0.9.1
    ## [13] pillar_1.6.4            lattice_0.20-45         glue_1.6.1
    ## [16] beachmat_2.10.0         digest_0.6.29           XVector_0.34.0
    ## [19] colorspace_2.0-2        cowplot_1.1.1           htmltools_0.5.2
    ## [22] pkgconfig_2.0.3         zlibbioc_1.40.0         purrr_0.3.4
    ## [25] scales_1.1.1            RSpectra_0.16-0         ScaledMatrix_1.2.0
    ## [28] BiocParallel_1.28.3     tibble_3.1.6            farver_2.1.0
    ## [31] generics_0.1.1          ellipsis_0.3.2          withr_2.4.3
    ## [34] magrittr_2.0.1          crayon_1.4.2            evaluate_0.14
    ## [37] fansi_1.0.2             FNN_1.1.3               beeswarm_0.4.0
    ## [40] tools_4.1.2             lifecycle_1.0.1         stringr_1.4.0
    ## [43] munsell_0.5.0           irlba_2.3.5             compiler_4.1.2
    ## [46] jquerylib_0.1.4         rsvd_1.0.5              rlang_0.4.12
    ## [49] grid_4.1.2              RCurl_1.98-1.5          BiocNeighbors_1.12.0
    ## [52] labeling_0.4.2          bitops_1.0-7            rmarkdown_2.11
    ## [55] gtable_0.3.0            DBI_1.1.2               R6_2.5.1
    ## [58] gridExtra_2.3           knitr_1.37              dplyr_1.0.7
    ## [61] uwot_0.1.11             fastmap_1.1.0           utf8_1.2.2
    ## [64] stringi_1.7.6           ggbeeswarm_0.6.0        parallel_4.1.2
    ## [67] Rcpp_1.0.8              vctrs_0.3.8             tidyselect_1.1.1
    ## [70] xfun_0.29               sparseMatrixStats_1.6.0