scRNA-seq稀疏矩阵图解,格式转换的核心

前言

  为什么单细胞数据要使用稀疏性矩阵?用一个实例来说明稀疏矩阵的效率用多好,且看下面的代码:

# 普通矩阵
dim(mat)
[1] 18104 23130

str(mat)
 num [1:18104, 1:23130] 0 0 0 0 0 0 0 0 0 0 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:18104(1d)] "FO538757.1" "AP006222.1" "AL732372.2" "AL669831.5" ...
  ..$ : chr [1:23130(1d)] "ccRCC_1_rep1@AAATGCCGTCTCAACA-1" "ccRCC_1_rep1@AACACGTAGATCACGG-1" "ccRCC_1_rep1@AACTTTCAGTCATCCA-1" "ccRCC_1_rep1@ACACCCTGTCGACTAT-1" ...

mat[19:26,58:59]
           ccRCC_1_rep2@CAGATCAAGCAATATG-1 ccRCC_1_rep2@CAGCGACCACAGACTT-1
TNFRSF4                          1.2251890                         0.00000
SDF4                             0.9567439                         0.00000
B3GALT6                          0.0000000                         0.00000
C1QTNF12                         0.0000000                         0.00000
AL162741.1                       0.0000000                         0.00000
UBE2J2                           0.5886769                         0.00000
SCNN1D                           0.0000000                         0.00000
ACAP3                            0.0000000                         1.58509

# 稀疏矩阵
dim(smat)
[1] 18104 23130

str(smat)
Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  ..@ i       : int [1:58298590] 19 31 37 58 81 98 109 237 248 264 ...
  ..@ p       : int [1:23131] 0 968 2263 3364 3919 4920 6049 7129 8248 9509 ...
  ..@ Dim     : int [1:2] 18104 23130
  ..@ Dimnames:List of 2
  .. ..$ : chr [1:18104(1d)] "FO538757.1" "AP006222.1" "AL732372.2" "AL669831.5" ...
  .. ..$ : chr [1:23130(1d)] "ccRCC_1_rep1@AAATGCCGTCTCAACA-1" "ccRCC_1_rep1@AACACGTAGATCACGG-1" "ccRCC_1_rep1@AACTTTCAGTCATCCA-1" "ccRCC_1_rep1@ACACCCTGTCGACTAT-1" ...
  ..@ x       : num [1:58298590] 1.54 1.54 1.54 1.54 2.13 ...
  ..@ factors : list()

smat[19:26,58:59]
8 x 2 sparse Matrix of class "dgCMatrix"
           ccRCC_1_rep2@CAGATCAAGCAATATG-1 ccRCC_1_rep2@CAGCGACCACAGACTT-1
TNFRSF4                          1.2251890                         .
SDF4                             0.9567439                         .
B3GALT6                          .                                 .
C1QTNF12                         .                                 .
AL162741.1                       .                                 .
UBE2J2                           0.5886769                         .
SCNN1D                           .                                 .
ACAP3                            .                                 1.58509

# 占用内存大小
print(object.size(mat),units='auto')
3.1 Gb
print(object.size(smat),units='auto')
670 Mb

  上面matsmat分别是普通矩阵和稀疏矩阵dgCMatrix,都是18104 x 23130的维度,使用起来也察觉不到异样。回到一开头的问题,那为什么单细胞要使用稀疏矩阵呢?答案就在上面的最后几行代码,可以看到存储同样大小的数据,稀疏矩阵占用的内存只约为普通矩阵的1/5,这效率就算是用笔记本,同时处理的样本数高低也能增加几个。
  与普通矩阵相比,稀疏矩阵的用法及视角效果虽然看上去没有什么区别,但是内部的存储方式却差别很大,以达到只保留非零值的目的,这也是效率差异的由来。虽说条条大路通罗马,但路与路不同,选择好走的路会让人事半功倍,效率大大的有。
  既然已经知道稀疏矩阵好使用,那自然会使用起来,可在此之前肯定要先创建矩阵,那么如何去创建一个稀疏矩阵呢?
  稀疏矩阵根据存储形式可以分成三种,如在R里面分别对应为dgTMatrixdgCMatrixdgRMatrix。其中,dgTMatrix原理最为简单直接,后面两个原理稍微拐了一个弯分别按照行和列压缩。

dgTMatrix

  该格式将普通矩阵的数据巧妙地一分为三(row,col,data),分别用来存储非零数值以及数值在矩阵中对应的行列索引。存储非零数值是所有稀疏矩阵格式的中心思想,删除大量的零值可不就提升了效率嘛。
  通常来说,单细胞数据中每个细胞里面表达的基因也就2000左右,按照编码基因20000个来算,也就1/10的基因有表达值,由此可见,稀疏矩阵和单细胞数据有点郎才女貌的感觉了,般配!

  从上面的示意图可以看出,每个数值在对应位置上都有自己的行列索引,所以这种方式存储的非零值顺序对结果没有影响

dgCMatrix

  dgCMatrix是一种按列压缩的稀疏矩阵格式,这应该是在dgTMatrix基础上演化而来。最大的区别在于非零值的存储是有顺序的,按列压缩的意思就是依据列的顺序,将第一列到最后一列的非零值依次存储起来,详见下面的示意图。

  这种方式也是将矩阵拆分成三部分,分别为indptrindicesdata,依据这三个相互配合提供的信息来定位非零值在矩阵的位置。

  • indptr:指示每列非零值的个数
  • indices:指示每个元素的行索引
  • data:元素的真实值

  indicesdata比较好理解,indptr并没有直接提供每列非零值的个数,而是拐了一个弯,用连续两个数的差(后面减去前面的差)做为个数。如上面的示意图,共5列,indptr却有6个值,按顺序分别求差就会得到1、0、3、2、1,分别指示1-5列非零值的个数。
  得到每列非零值的个数后,接着从indices里面顺序寻找相应个数的索引,这些索引分别指示每个元素在矩阵的行位置。如第一列有1个非零值,然后分别从indicesdata里面找这个非零值的行索引和真实值,最终确定矩阵中第一列第一行的值为8。

dgRMatrix

  该格式与dgCMatrix原理一致,只不过是将稀疏矩阵按行压缩,也就是说indptr变成了指示每行非零值的个数,而indices变成了指示每个元素的列索引。

结束语

  格式转换的核心就是知道数据的内部结构。为什么内部原理很重要?因为眼睛会欺骗大脑,这也用一句名言来解释:眼见不一定为实。就像稀疏矩阵,看起来与普通矩阵无异,但内部存储结构却截然不同。所以,了解了内部原理,不论给予什么样的数据,咱们都可以拼凑出想要的数据格式,比如h5loom等格式读取为矩阵,感兴趣可以阅读具体之前的帖子[scRNAseq:h5文件转化为matrix表达矩阵]、[一网打尽scRNA矩阵格式读取和转化(h5 h5ad loom)]。

参考

  做为网络拾遗 (个人觉得用这个词,来形容网络冲浪时吸收到的内容,挺贴切的,毕竟不是有句说得好:读书人的事,能叫偷吗!) 的收获,还是要感谢博主的分享:https://blog.csdn.net/jeffery0207/article/details/122507934


往期回顾

venn | 多样本间peak重叠韦恩图的解决方案
【海外招聘】美国威尔康奈尔医学院招聘博士后
R编程技巧 | 学习高手实现函数多功能化的两种方法
linux | while + read 实现本地 or 集群批处理,实用且优雅
pyscenic | 单细胞转录因子分析,原理图文详解

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 218,451评论 6 506
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 93,172评论 3 394
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 164,782评论 0 354
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 58,709评论 1 294
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 67,733评论 6 392
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 51,578评论 1 305
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 40,320评论 3 418
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 39,241评论 0 276
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,686评论 1 314
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,878评论 3 336
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,992评论 1 348
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,715评论 5 346
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 41,336评论 3 330
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,912评论 0 22
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 33,040评论 1 270
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 48,173评论 3 370
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,947评论 2 355

推荐阅读更多精彩内容