# 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 ```