单细胞交响乐
  • 前言:我与《单细胞交响乐》的缘分
  • 1 准备篇:背景知识
    • 1.1 数据结构
    • 1.2 总览 | 从实验到分析
  • 2 积累篇:文献阅读
    • 2.1.1 综述 | 2019-单细胞转录组分析最佳思路
    • 2.1.2 综述 | 2018-单细胞捕获平台
    • 2.1.3 综述 | 2017-scRNA中的细胞聚类分群
    • 2.1.4 综述 | scRNA已经开发出超过1000款工具了,你用过几种?
    • 2.1.5 综述 | 2021-单细胞测序的微流控技术应用
    • 2.2.1 研究 | 2018-单细胞转录组探索癌症免疫治疗获得性抗性机理
    • 2.2.2 研究 | 2018-人类结直肠癌单细胞多组学分析
    • 2.2.3 研究 | 2020-单细胞分析揭示葡萄膜黑色素瘤新的进化复杂性
    • 2.2.4 研究 | 2020-COVID-19病人支气管免疫细胞单细胞测序分析
    • 2.2.5 研究 | 2020-原汁原味读--单细胞肿瘤免疫图谱
    • 2.2.6 研究 | 2021-多发性骨髓瘤发展过程中肿瘤和免疫细胞的共同进化
    • 2.2.7 研究 | 2021-多个组织的成纤维细胞图谱
    • 2.2.8 研究 | 2021-多组学分析肺结核队列的记忆T细胞状态
    • 2.2.9 研究 | 2021-CancerSCEM: 人类癌症单细胞表达图谱数据库
    • 2.2.10 研究| 2021-单细胞转录组分析COVID-19重症患者肺泡巨噬细胞亚型
    • 2.2.11 研究 |2021-单细胞转录组揭示肺腺癌特有的肿瘤微环境
    • 2.2.12 研究 | 2021-单细胞转录组揭示乳头状甲状腺癌起始与发展
    • 2.2.13 研究 | 2021-解析食管鳞癌化疗病人的单细胞转录组
    • 2.2.14 研究 | 2021-单细胞水平看骨髓瘤的细胞状态和基因调控
    • 2.3.1 算法|2020-BatchBench比较scRNA批次矫正方法
    • 2.3.2 算法 | 2021-scPhere——用地球仪来展示降维结果
    • 2.3.3 算法 | 2021-单细胞差异分析方法评测
    • 2.3.4 算法 | 2021-细胞分群新方法——CNA(co-varying neighborhood analysis)
    • 2.3.5 工具 | 2018-iSEE:单细胞数据可视化辅助网页工具
    • 2.3.6 工具 | 2021-MACA: 一款自动注释细胞类型的工具
    • 2.3.7 工具 | 2021-一个很有想法的工具——Ikarus,想要在单细胞水平直接鉴定肿瘤细胞
  • 3 流程篇:分析框架
    • 3.1 质控
    • 3.2 归一化
    • 3.3 挑选表达量高变化基因
    • 3.4 降维
    • 3.5 聚类
    • 3.6 Marker/标记基因检测
    • 3.7 细胞类型注释
    • 3.8 批次效应处理
    • 3.9 多样本间差异分析
    • 3.10 检测Doublet
    • 3.11 细胞周期推断
    • 3.12 细胞轨迹推断
    • 3.13 与蛋白丰度信息结合
    • 3.14 处理大型数据
    • 3.15 不同R包数据的相互转换
  • 4 实战篇:活学活用
    • 4.1 实战一 | Smart-seq2 | 小鼠骨髓
    • 4.2 实战二 | STRT-Seq | 小鼠大脑
    • 4.3 实战三 | 10X | 未过滤的PBMC
    • 4.4 实战四 | 10X | 过滤后的PBMC
    • 4.5 实战五 | CEL-seq2 | 人胰腺细胞
    • 4.6 实战六 | CEL-seq | 人胰腺细胞
    • 4.7 实战七 | SMARTer | 人胰腺细胞
    • 4.8 实战八 | Smart-seq2 | 人胰腺细胞
    • 4.9 实战九 | 不同技术数据整合 | 人胰腺细胞
    • 4.10 实战十 | CEL-seq | 小鼠造血干细胞
    • 4.11 实战十一 | Smart-seq2 | 小鼠造血干细胞
    • 4.12 实战十二 | 10X | 小鼠嵌合体胚胎
    • 4.13 实战十三 | 10X | 小鼠乳腺上皮细胞
    • 4.14 | 实战十四 | 10X | HCA计划的38万骨髓细胞
  • 5 补充篇:开拓思路
    • 5.1 10X Genomics概述
      • 5.1.1 10X Genomics 问题集锦
    • 5.2 CellRanger篇
      • 5.2.1 CellRanger实战(一)数据下载
      • 5.2.2 CellRanger实战(二) 使用前注意事项
      • 5.2.3 CellRanger实战(三) 使用初探
      • 5.2.4 CellRanger实战(四)流程概览
      • 5.2.5 CellRanger实战(五) 理解count输出的结果
    • 5.3 Seurat的使用
      • 5.3.1 Seurat V3 | 实战之2700 PBMCs分析
      • 5.3.2 Seurat V3 | 如何改造Seurat包的DoHeatmap函数?
      • 5.3.3 scRNA的3大R包对比
      • 5.3.4 Seurat两种数据比较:integrated vs RNA assay
      • 5.3.5 seurat 的几种findmaker比较
    • 5.4 Monocle的使用
      • 5.4.1 Monocle V3实战
    • 5.5 多个数据集的整合
      • 5.5.1 使用Seurat的merge功能进行整合
      • 5.5.2 如何使用sctransform去除批次效应
由 GitBook 提供支持
在本页
  • 1 前言
  • 2 快速估算
  • 2.1 近似而非精确的近邻搜索
  • 2.2 奇异值分解
  • 3 并行计算
  • 3.1 为什么?
  • 3.2 怎么做?
  • 注意
  • 4 可能会遇到内存不足

这有帮助吗?

  1. 3 流程篇:分析框架

3.14 处理大型数据

刘小泽写于2020.7.18

1 前言

scRNA-seq技术越来越成熟,测的细胞数也在日益增长,越来越多的研究结果发表在GEO数据库中,而且还有大型单细胞项目的开展(例如人类图谱计划HCA)。因此分析方法需要考虑到逐渐庞大的数据集,这次就来看看怎么做可以帮助我们获得更快的处理速度和分析效率。

2 快速估算

2.1 近似而非精确的近邻搜索

在高维空间中对近邻细胞进行判断基本上是必备流程,像buildSNNGraph(), doubletCells() 都是需要经历这一步。默认是利用KNN(k-nearest neighbours)算法找到更准确的近邻数量,而牺牲速度。但是大型数据更需要的是速度,例如这个BiocNeighbors R包就可以通过BNPARAM=轻松转换近邻搜索方法

以pbmc数据为例,这个数据之前应该分享过

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

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)

这里因为数据量不大,并看不出速度上的优势,只是为了证明二者在准确度上相差无几

2.2 奇异值分解

术语叫做:singular value decomposition (SVD),在 denoisePCA(), fastMNN(), doubletCells()都有应用。默认使用base::svd() 也是进行了精确的计算,同样不适合大型数据集。好在 irlba和rsvd 两个R包都提供了近似计算方法,具体的参数配置和AnnoyParam 一样也是做成了BiocSingular包的一个函数

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的分析中

不同的硬件和操作系统,选择的并行方法也不同

library(BiocParallel)
registered() #可以看看支持哪些类型的加速,最顶上的是默认的

一般包括:

  • SerialParam:全平台支持

  • MulticoreParam:适用于Unix和Mac。在windows上它等同于SerialParam

  • SnowParam:全平台支持

  • BatchtoolsParam:集群

  • DoparParam:全平台支持

如果要更改

# 本来的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)

可以这么使用

# 不同的方法
dec.pbmc.mc <- modelGeneVar(sce.pbmc, BPPARAM=MulticoreParam(2))
dec.pbmc.snow <- modelGeneVar(sce.pbmc, BPPARAM=SnowParam(5))

在SLURM HPC中,可以使用BatchtoolsParam

# 设置每个任务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,就是用硬盘空间来存储数据,必要时再调用一部分数据到内存,虽然这一来一回很影响处理速度,但毕竟可以用。

使用这个HDF5ArrayR包可以做到(类似的还有 bigmemory, matter),它会将底层数据做成HDF5格式

例如,从130万个脑细胞数据中选取了2万个

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的矩阵

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

看一下这个大小

object.size(counts(sce.brain))
## 2328 bytes

但实际上这个底层数据的大小是

file.info(path(counts(sce.brain)))$size
## [1] 76264332
上一页3.13 与蛋白丰度信息结合下一页3.15 不同R包数据的相互转换

最后更新于4年前

这有帮助吗?

参考:

参考:

参考:

https://cosx.org/2016/09/r-and-parallel-computing/
https://bioconductor.org/packages/3.11/bioc/vignettes/BiocParallel/inst/doc/Introduction_To_BiocParallel.pdf
https://bioconductor.org/packages/3.11/BiocParallel/vignettes/BiocParallel_BatchtoolsParam.pdf