• 主页
  • 课程

    关于课程

    • 课程归档
    • 成为一名讲师
    • 讲师信息
    同等学历教学

    同等学历教学

    免费
    阅读更多
  • 特色
    • 展示
    • 关于我们
    • 问答
  • 事件
  • 个性化
  • 博客
  • 联系
  • 站点资源
    有任何问题吗?
    (00) 123 456 789
    weinfoadmin@weinformatics.cn
    注册登录
    恒诺新知
    • 主页
    • 课程

      关于课程

      • 课程归档
      • 成为一名讲师
      • 讲师信息
      同等学历教学

      同等学历教学

      免费
      阅读更多
    • 特色
      • 展示
      • 关于我们
      • 问答
    • 事件
    • 个性化
    • 博客
    • 联系
    • 站点资源

      生信星球

      • 首页
      • 博客
      • 生信星球
      • 单细胞交响乐-处理大型数据

      单细胞交响乐-处理大型数据

      • 发布者 weinfoeditor
      • 分类 生信星球
      • 日期 2020年7月24日
      测试开头

      单细胞交响乐-处理大型数据单细胞交响乐-处理大型数据 今天是生信星球陪你的第685天单细胞交响乐-处理大型数据


         大神一句话,菜鸟跑半年。我不是大神,但我可以缩短你走弯路的半年~

         就像歌儿唱的那样,如果你不知道该往哪儿走,就留在这学点生信好不好~

         这里有豆豆和花花的学习历程,从新手到进阶,生信路上有你有我!

      豆豆写于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 为什么?

      参考:https://cosx.org/2016/09/r-and-parallel-computing/

      R 采用的是内存计算模式(In-Memory),被处理的数据需要预取到主存(RAM)中。其优点是计算效率高、速度快,但缺点是这样一来能处理的问题规模就非常有限(小于 RAM 的大小)。另一方面,R 的核心(R core)是一个单线程的程序,在多核处理器上,R 无法有效地利用所有的计算内核。即使机器性能很强大,有32个核心,但它也只能使用1/32的计算能力,浪费了31/32。

      3.2 怎么做?

      使用 BiocParallel 包,可以将并行运算覆盖到基于Bioconductor的分析中

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

      参考:
      https://bioconductor.org/packages/3.11/bioc/vignettes/BiocParallel/inst/doc/Introduction_To_BiocParallel.pdf

      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

      参考:
      https://bioconductor.org/packages/3.11/BiocParallel/vignettes/BiocParallel_BatchtoolsParam.pdf

      # 设置每个任务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

      点击底部的“阅读原文”,获得更好的阅读体验哦😻

      初学生信,很荣幸带你迈出第一步

      🤓生信星球 🌎~ 一个不拽术语、通俗易懂的生信知识平台

      测试结尾

      请关注“恒诺新知”微信公众号,感谢“R语言“,”数据那些事儿“,”老俊俊的生信笔记“,”冷🈚️思“,“珞珈R”,“生信星球”的支持!

      • 分享:
      作者头像
      weinfoeditor

      上一篇文章

      找不到某些GEO数据的表达矩阵肿么办
      2020年7月24日

      下一篇文章

      生信爆款入门-第7期(线上直播4周,马拉松式陪伴,带你入门)
      2020年7月26日

      你可能也喜欢

      8-1651673488
      生信零基础入门学习小组长期报名中(2022仍继续)
      7 4月, 2022
      2-1651673738
      简化版的ROC曲线
      21 2月, 2022
      8-1651674718
      支持向量机模型
      19 11月, 2021

      搜索

      分类

      • R语言
      • TCGA数据挖掘
      • 单细胞RNA-seq测序
      • 在线会议直播预告与回放
      • 数据分析那些事儿分类
      • 未分类
      • 生信星球
      • 老俊俊的生信笔记

      投稿培训

      免费

      alphafold2培训

      免费

      群晖配置培训

      免费

      最新博文

      白介素-17受体信号的自主激活,维持炎症并促进疾病进展
      048月2023
      MCT4依赖的乳酸分泌抑制LKB1缺陷肺腺癌的抗肿瘤免疫
      187月2023
      Nature | 单细胞技术揭示衰老细胞与肌肉再生
      301月2023
      logo-eduma-the-best-lms-wordpress-theme

      (00) 123 456 789

      weinfoadmin@weinformatics.cn

      恒诺新知

      • 关于我们
      • 博客
      • 联系
      • 成为一名讲师

      链接

      • 课程
      • 事件
      • 展示
      • 问答

      支持

      • 文档
      • 论坛
      • 语言包
      • 发行状态

      推荐

      • iHub汉语代码托管
      • iLAB耗材管理
      • WooCommerce
      • 丁香园论坛

      weinformatics 即 恒诺新知。ICP备案号:粤ICP备19129767号

      • 关于我们
      • 博客
      • 联系
      • 成为一名讲师

      要成为一名讲师吗?

      加入数以千计的演讲者获得100%课时费!

      现在开始

      用你的站点账户登录

      忘记密码?

      还不是会员? 现在注册

      注册新帐户

      已经拥有注册账户? 现在登录

      close
      会员购买 你还没有登录,请先登录
      • ¥99 VIP-1个月
      • ¥199 VIP-半年
      • ¥299 VIP-1年
      在线支付 激活码

      立即支付
      支付宝
      微信支付
      请使用 支付宝 或 微信 扫码支付
      登录
      注册|忘记密码?