前言
为什么单细胞数据要使用稀疏性矩阵？用一个实例来说明稀疏矩阵的效率用多好，且看下面的代码：
# 普通矩阵
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
上面mat和smat分别是普通矩阵和稀疏矩阵dgCMatrix，都是18104 x 23130的维度，使用起来也察觉不到异样。回到一开头的问题，那为什么单细胞要使用稀疏矩阵呢？答案就在上面的最后几行代码，可以看到存储同样大小的数据，稀疏矩阵占用的内存只约为普通矩阵的1/5，这效率就算是用笔记本，同时处理的样本数高低也能增加几个。
与普通矩阵相比，稀疏矩阵的用法及视角效果虽然看上去没有什么区别，但是内部的存储方式却差别很大，以达到只保留非零值的目的，这也是效率差异的由来。虽说条条大路通罗马，但路与路不同，选择好走的路会让人事半功倍，效率大大的有。
既然已经知道稀疏矩阵好使用，那自然会使用起来，可在此之前肯定要先创建矩阵，那么如何去创建一个稀疏矩阵呢？
稀疏矩阵根据存储形式可以分成三种，如在R里面分别对应为dgTMatrix、dgCMatrix、dgRMatrix。其中，dgTMatrix原理最为简单直接，后面两个原理稍微拐了一个弯分别按照行和列压缩。
dgTMatrix
该格式将普通矩阵的数据巧妙地一分为三(row，col，data)，分别用来存储非零数值以及数值在矩阵中对应的行列索引。存储非零数值是所有稀疏矩阵格式的中心思想，删除大量的零值可不就提升了效率嘛。
通常来说，单细胞数据中每个细胞里面表达的基因也就2000左右，按照编码基因20000个来算，也就1/10的基因有表达值，由此可见，稀疏矩阵和单细胞数据有点郎才女貌的感觉了，般配！
从上面的示意图可以看出，每个数值在对应位置上都有自己的行列索引，所以这种方式存储的非零值顺序对结果没有影响。
dgCMatrix
dgCMatrix是一种按列压缩的稀疏矩阵格式，这应该是在dgTMatrix基础上演化而来。最大的区别在于非零值的存储是有顺序的，按列压缩的意思就是依据列的顺序，将第一列到最后一列的非零值依次存储起来，详见下面的示意图。
这种方式也是将矩阵拆分成三部分，分别为indptr、indices、data，依据这三个相互配合提供的信息来定位非零值在矩阵的位置。
-
indptr：指示每列非零值的个数 -
indices：指示每个元素的行索引 -
data：元素的真实值
indices和data比较好理解，indptr并没有直接提供每列非零值的个数，而是拐了一个弯，用连续两个数的差(后面减去前面的差)做为个数。如上面的示意图，共5列，indptr却有6个值，按顺序分别求差就会得到1、0、3、2、1，分别指示1-5列非零值的个数。
得到每列非零值的个数后，接着从indices里面顺序寻找相应个数的索引，这些索引分别指示每个元素在矩阵的行位置。如第一列有1个非零值，然后分别从indices和data里面找这个非零值的行索引和真实值，最终确定矩阵中第一列第一行的值为8。
dgRMatrix
该格式与dgCMatrix原理一致，只不过是将稀疏矩阵按行压缩，也就是说indptr变成了指示每行非零值的个数，而indices变成了指示每个元素的列索引。
结束语
格式转换的核心就是知道数据的内部结构。为什么内部原理很重要？因为眼睛会欺骗大脑，这也用一句名言来解释：眼见不一定为实。就像稀疏矩阵，看起来与普通矩阵无异，但内部存储结构却截然不同。所以，了解了内部原理，不论给予什么样的数据，咱们都可以拼凑出想要的数据格式，比如h5、loom等格式读取为矩阵，感兴趣可以阅读具体之前的帖子[scRNAseq：h5文件转化为matrix表达矩阵]、[一网打尽scRNA矩阵格式读取和转化(h5 h5ad loom)]。
参考
做为网络拾遗 (个人觉得用这个词，来形容网络冲浪时吸收到的内容，挺贴切的，毕竟不是有句说得好：读书人的事，能叫偷吗！) 的收获，还是要感谢博主的分享：https://blog.csdn.net/jeffery0207/article/details/122507934。
往期回顾
venn | 多样本间peak重叠韦恩图的解决方案
【海外招聘】美国威尔康奈尔医学院招聘博士后
R编程技巧 | 学习高手实现函数多功能化的两种方法
linux | while + read 实现本地 or 集群批处理，实用且优雅
pyscenic | 单细胞转录因子分析，原理图文详解
网友评论