# 3.14 处理大型数据 ## 1 前言 scRNA-seq技术越来越成熟，测的细胞数也在日益增长，越来越多的研究结果发表在GEO数据库中，而且还有大型单细胞项目的开展（例如人类图谱计划HCA）。因此分析方法需要考虑到逐渐庞大的数据集，这次就来看看怎么做可以帮助我们获得更快的处理速度和分析效率。 ## 2 快速估算 ### **2.1 近似而非精确的近邻搜索** 在高维空间中对近邻细胞进行判断基本上是必备流程，像`buildSNNGraph()`, `doubletCells()` 都是需要经历这一步。默认是利用KNN（k-nearest neighbours）算法找到更准确的近邻数量，而牺牲速度。但是大型数据更需要的是速度，例如这个`BiocNeighbors` R包就可以通过`BNPARAM=`轻松转换近邻搜索方法 以pbmc数据为例，这个数据之前应该分享过 ```r load('clustered.sce.pbmc.RData') sce.pbmc ## class: SingleCellExperiment ## dim: 33694 3922 ## metadata(1): Samples ## assays(2): counts logcounts ## rownames(33694): RP11-34P13.3 FAM138A ... AC213203.1 FAM231B ## rowData names(2): ID Symbol ## colnames(3922): AAACCTGAGAAGGCCT-1 AAACCTGAGACAGACC-1 ... ## TTTGTCACAGGTCCAC-1 TTTGTCATCCCAAGAT-1 ## colData names(4): Sample Barcode sizeFactor label ## reducedDimNames(3): PCA TSNE UMAP ## altExpNames(0): ``` 看到这里的数据是经过了PCA、t-SNE、UMAP降维操作的，并已经做好了分群，使用的也是默认的精确KNN搜索。 **下面使用近似搜索：** `AnnoyParam`的解释是：A class to hold parameters for the Annoy algorithm for approximate nearest neighbor identification. 也就是它包装了近似搜索的一套参数，并传递给`buildSNNGraph` ```r library(scran) library(BiocNeighbors) snn.gr <- buildSNNGraph(sce.pbmc, BNPARAM=AnnoyParam(), use.dimred="PCA") clusters <- igraph::cluster_walktrap(snn.gr) table(Exact=colLabels(sce.pbmc), Approx=clusters$membership) ``` ![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-18-005159.png) 这里因为数据量不大，并看不出速度上的优势，只是为了证明二者在准确度上相差无几 ### **2.2 奇异值分解** 术语叫做：singular value decomposition (SVD)，在 `denoisePCA()`, `fastMNN()`, `doubletCells()`都有应用。默认使用`base::svd()` 也是进行了精确的计算，同样不适合大型数据集。好在 `irlba`和`rsvd` 两个R包都提供了近似计算方法，具体的参数配置和`AnnoyParam` 一样也是做成了`BiocSingular`包的一个函数 ```r library(scater) library(BiocSingular) # 方法一：randomized SVD (RSVD) set.seed(101000) r.out <- runPCA(sce.pbmc, ncomponents=20, BSPARAM=RandomParam()) str(reducedDim(r.out)) ## num [1:3922, 1:20] 15.3 13.41 -8.46 -7.86 6.38 ... ## - attr(*, "dimnames")=List of 2 ## ..$ : chr [1:3922] "AAACCTGAGAAGGCCT-1" "AAACCTGAGACAGACC-1" "AAACCTGAGGCATGGT-1" "AAACCTGCAAGGTTCT-1" ... ## ..$ : chr [1:20] "PC1" "PC2" "PC3" "PC4" ... ## - attr(*, "percentVar")= num [1:20] 20.26 10.02 5.36 2.19 1.41 ... ## - attr(*, "rotation")= num [1:500, 1:20] 0.2015 0.182 0.1764 0.1067 0.0649 ... ## ..- attr(*, "dimnames")=List of 2 ## .. ..$ : chr [1:500] "LYZ" "S100A9" "S100A8" "HLA-DRA" ... ## .. ..$ : chr [1:20] "PC1" "PC2" "PC3" "PC4" ... # 方法二：IRLBA set.seed(101001) i.out <- runPCA(sce.pbmc, ncomponents=20, BSPARAM=IrlbaParam()) str(reducedDim(i.out)) ## num [1:3922, 1:20] 15.3 13.41 -8.46 -7.86 6.38 ... ## - attr(*, "dimnames")=List of 2 ## ..$ : chr [1:3922] "AAACCTGAGAAGGCCT-1" "AAACCTGAGACAGACC-1" "AAACCTGAGGCATGGT-1" "AAACCTGCAAGGTTCT-1" ... ## ..$ : chr [1:20] "PC1" "PC2" "PC3" "PC4" ... ## - attr(*, "percentVar")= num [1:20] 20.26 10.02 5.36 2.19 1.41 ... ## - attr(*, "rotation")= num [1:500, 1:20] 0.2015 0.182 0.1764 0.1067 0.0649 ... ## ..- attr(*, "dimnames")=List of 2 ## .. ..$ : chr [1:500] "LYZ" "S100A9" "S100A8" "HLA-DRA" ... ## .. ..$ : chr [1:20] "PC1" "PC2" "PC3" "PC4" ... ``` 这两种方法的速度都比准确的SVD方法快，丢失的准确度也微乎其微，也因此被scran和scater包设为了默认的参数，也因此需要设置随机种子（毕竟都是估计的方法，每次运行结果都不一致）。当然，**IRLBA相比RSVD更准确，但速度不如RSVD。** ## 3 并行计算 ### **3.1 为什么？** > 参考： R 采用的是内存计算模式（In-Memory），被处理的数据需要预取到主存（RAM）中。其优点是计算效率高、速度快，但缺点是这样一来能处理的问题规模就非常有限（小于 RAM 的大小）。另一方面，R 的核心（R core）是一个单线程的程序，在多核处理器上，R 无法有效地利用所有的计算内核。即使机器性能很强大，有32个核心，但它也只能使用1/32的计算能力，浪费了31/32。 ### **3.2 怎么做？** **使用BiocParallel包，可以将并行运算覆盖到基于Bioconductor的分析中** 不同的硬件和操作系统，选择的并行方法也不同 > 参考： ```r library(BiocParallel) registered() #可以看看支持哪些类型的加速，最顶上的是默认的 ``` 一般包括： * SerialParam：全平台支持 * MulticoreParam：适用于Unix和Mac。在windows上它等同于SerialParam * SnowParam：全平台支持 * BatchtoolsParam：集群 * DoparParam：全平台支持 如果要更改 ```r # 本来的default是 MulticoreParam default <- registered() # 现在改成BatchtoolsParam register(BatchtoolsParam(workers = 10), default = TRUE) names(registered()) ## [1] "BatchtoolsParam" "MulticoreParam" "SnowParam" "SerialParam" # 要再恢复原来的设置 for (param in rev(default)) register(param) ``` 可以这么使用 ```r # 不同的方法 dec.pbmc.mc <- modelGeneVar(sce.pbmc, BPPARAM=MulticoreParam(2)) dec.pbmc.snow <- modelGeneVar(sce.pbmc, BPPARAM=SnowParam(5)) ``` 在SLURM HPC中，可以使用`BatchtoolsParam` > 参考： ```r # 设置每个任务2小时、8G内存、1CPU，总共10个任务 bpp <- BatchtoolsParam(10, cluster="slurm", resources=list(walltime=7200, memory=8000, ncpus=1)) ``` ### **注意** * 这个并行计算只是加速了CPU的计算速度，但如果任务受限于内存或硬盘的读入读出，它依然是没办法加速的； * R实现多线程还是很麻烦的，它的操作逻辑是：先设置一个或多个session（对话窗口）；然后加载相关的包；session之间进行数据传递。可能最后还不如单线程运行效果好 ## 4 可能会遇到内存不足 想一下，我们处理的核心是不是基于表达矩阵？既然是核心，就需要完全载入内存中，然后才能实现后面的顺利读取、处理。现在单细胞的表达矩阵有两种形式：稀疏矩阵`dgCMatrix`或者普通矩阵`matrix` 。如果我们的数据是10X产生的130万个脑细胞数据，如果是以普通矩阵读入，就需要消耗100G内存，即使使用稀疏矩阵，也需要消耗30G内存左右。因此没有很好的服务器基本上干不了这个事。 **用处理速度换取内存不足** 如果真的面临无法增加内存的困境，还有一个plan B，就是用硬盘空间来存储数据，必要时再调用一部分数据到内存，虽然这一来一回很影响处理速度，但毕竟可以用。 使用这个`HDF5Array`R包可以做到（类似的还有 `bigmemory`, `matter`），它会将底层数据做成HDF5格式 例如，从130万个脑细胞数据中选取了2万个 ```r library(TENxBrainData) sce.brain <- TENxBrainData20k() sce.brain ## class: SingleCellExperiment ## dim: 27998 20000 ## metadata(0): ## assays(1): counts ## rownames: NULL ## rowData names(2): Ensembl Symbol ## colnames: NULL ## colData names(4): Barcode Sequence Library Mouse ## reducedDimNames(0): ## altExpNames(0): ``` 看一下这个表达矩阵，显然是一个HDF5的矩阵 ```r counts(sce.brain) ## <27998 x 20000> matrix of class HDF5Matrix and type "integer": ## [,1] [,2] [,3] [,4] ... [,19997] [,19998] [,19999] ## [1,] 0 0 0 0 . 0 0 0 ## [2,] 0 0 0 0 . 0 0 0 ## [3,] 0 0 0 0 . 0 0 0 ## [4,] 0 0 0 0 . 0 0 0 ## [5,] 0 0 0 0 . 0 0 0 ## ... . . . . . . . . ## [27994,] 0 0 0 0 . 0 0 0 ## [27995,] 0 0 0 1 . 0 2 0 ## [27996,] 0 0 0 0 . 0 1 0 ## [27997,] 0 0 0 0 . 0 0 0 ## [27998,] 0 0 0 0 . 0 0 0 ## [,20000] ## [1,] 0 ## [2,] 0 ## [3,] 0 ## [4,] 0 ## [5,] 0 ## ... . ## [27994,] 0 ## [27995,] 0 ## [27996,] 0 ## [27997,] 0 ## [27998,] 0 ``` 看一下这个大小 ```r object.size(counts(sce.brain)) ## 2328 bytes ``` 但实际上这个底层数据的大小是 ```r file.info(path(counts(sce.brain)))$size ## [1] 76264332 ``` --- # Agent Instructions: Querying This Documentation If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question. Perform an HTTP GET request on the current page URL with the `ask` query parameter: ``` GET https://jieandze1314.osca.top/03/03-14.md?ask= ``` The question should be specific, self-contained, and written in natural language. The response will contain a direct answer to the question and relevant excerpts and sources from the documentation. Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.