稀疏矩阵几乎产生于所有的大型科学工程计算领域，记录样本特征值的稠密矩阵中很多记录值都是0，由于0不携带信息，因此耗费空间存储0元素是很浪费资源的行为。而且很多计算只对非零元素进行操作，将特征矩阵构建成稀疏矩阵，可以很容易的索引到非零元素，所以基于稀疏矩阵的数据运算，可是实现更低的资源占用和更快的计算速度。
在单细胞领域，稀疏矩阵对于处理 scRNA-seq 表达谱数据是非常必要的，构建分析对象的时候它们消耗更低的内存。本文重点介绍 R语言 稀疏矩阵对象格式，稀疏矩阵的重构， 稀疏矩阵与稠密矩阵的相互转换。
三种类型的稀疏矩阵对象结构
首先，稀疏矩阵的一般格式主要有三种：
-
dgTMatrix ::: 创建未压缩的 三元组 稀疏矩阵 - 坐标格式的稀疏矩阵 - coo
coo格式的稀疏矩阵以三元组(i,j,x)的形式记录非零条目-
i## 非零元素的行号【0-based索引值】 -
j## 非零元素的列号【0-based索引值】 -
x## 按列向量方向排列的非零元素 值（从左到右）
-
-
dgCMatrix ::: 创建列压缩的 三元组 稀疏矩阵 - 压缩系数矩阵 - csc - 更快的列检索速度
csc格式的稀疏矩阵以三元组(i,p,x)的形式记录非零条目-
i## 非零元素的行号【0-based索引值】 -
p## 从0开始, 记录值是对每列非零元素的个数的累加（从左到右） -
x## 按列向量方向排列的非零元素 值 （从左到右）
-
-
dgRMatrix ::: 创建行压缩的 三元组 稀疏矩阵 - csr - 更快的行检索速度
-
p## 从0开始, 记录值是对每行非零元素的个数的累加（从上到下） -
j## 非零元素的列号【0-based索引值】 -
x## 按行向量方向排列的非零元素 值 （从上到下）
-
1、R 语言 Matrix包
Matrix 包为扩展基本矩阵数据类型的密集和稀疏矩阵提供了一组 S4 类和应用于这些类对象的各种函数和运算符方法。>Matrix package - RDocumentation
install.packages("Matrix") ; library(Matrix)
2、R 语言 稠密矩阵 与 稀疏矩阵 的相互转换
构建一个测试用稠密矩阵：
library(dplyr)
counts <- data.frame(cell.1=c(2,0,0,8), cell.2=c(4,4,0,9),cell.3=c(3,3,3,0),cell.4=c(3,3,3,3),row.names = paste0("gene.",seq(4))) %>% as.matrix()
> counts
cell.1 cell.2 cell.3 cell.4
gene.1 2 4 3 3
gene.2 0 4 3 3
gene.3 0 0 3 3
gene.4 8 9 0 3
> print( paste0("矩阵非零元素个数 ", sum( colSums(counts != 0)) ," 个") )
[1] "矩阵非零元素个数 12 个"
> is(counts, 'sparseMatrix') #判断为非 稀疏矩阵 对象
[1] FALSE
- 通过
as( )方法转换 稠密矩阵稀疏矩阵
sparse.gbm.T <- as(counts, "dgTMatrix") ### convert to coo_matrix
sparse.gbm.C <- as(counts, "dgCMatrix") ### convert to csc_matrix
sparse.gbm.R <- as(counts, "dgRMatrix") ### convert to csr_matrix
创建列压缩的稀疏矩阵(csc_matrix)还可以用方法sparse.gbm <- Matrix::Matrix(counts, sparse = T)，但是该方法在大型矩阵中较为低效。
三种稀疏矩阵对象在R语言中的结构str(sparse.gbm.*)
数据结构
关于压缩格式的稀疏矩阵对象中数据结构的 sparse.gbm.C@p 向量，该向量不是非零元素的行/列索引编号，它的值是对 数据压缩方向上 每个向量内非零元素个数的累加(如创建列压缩稀疏矩阵时，压缩方向为按列 向左到右 移动)。p向量的记录方式相比i/j记录非零元素索引号，有助于进一步减少稀疏矩阵对象存储大小。
p的计算
对于任意矩阵，都可以计算
p列向量的p c(0, cumsum(colSums(counts != 0)))行向量的p c(0, cumsum(rowSums(counts != 0)))通过
p 可以解析出 压缩方向的 真是索引编号(i/j)
#通过 p向量 解析出 列压缩稀疏矩阵的 向量j [记录了每个非零元素的列编号]
> rep(1:sparse.gbm.C@Dim[2], diff(sparse.gbm@p))
[1] 1 1 2 2 2 3 3 3 4 4 4 4
> sparse.gbm.C@i
[1] 0 3 0 1 3 0 1 2 0 1 2 3
> sparse.gbm.C@x
[1] 2 8 4 4 9 3 3 3 3 3 3 3
> sparse.gbm.C@j
Error: no slot of name "j" for this object of class "dgCMatrix"
#列压缩的稀疏矩阵用p向量记录了元素的列索引编号j，所以没有向量j
- 通过
as( )方法转换 稀疏矩阵
>as.matrix(sparse.gbm.T)
cell.1 cell.2 cell.3 cell.4
gene.1 2 4 3 3
gene.2 0 4 3 3
gene.3 0 0 3 3
gene.4 8 9 0 3
3、重构sparseMatrix
sparse.gbm.C.new <- Matrix::sparseMatrix(
i = sparse.gbm.C@i + 1 , # +1 转换 0-based编号 为R语言 1-based 编号方式
p = sparse.gbm.C@p ,
x = sparse.gbm.C@x,
dims = sparse.gbm.C@Dim #矩阵的维度
)
由于dgCMatrix对象的行/列索引值 i / j 是【0-based】的，而R中数据对象（如data.frame,vector）的索引值都是 【1-based】 的，所以在通过Matrix::sparseMatrix重构稀疏矩阵的时候，i/j的索引值要加1转换为【1-based】
除了按S4的方法从对象中取三元组信息，通过Matrix::summary()方法可以快速获取稀疏矩阵S4对象的（i,j,x)三元组，并用于重构 sparse_matrix。
> as.data.frame(summary(sparse.gbm.C))
i j x
1 1 1 2
2 4 1 8
3 1 2 4
4 2 2 4
5 4 2 9
6 1 3 3
7 2 3 3
8 3 3 3
9 1 4 3
10 2 4 3
11 3 4 3
12 4 4 3
4、稀疏矩阵的一些运算
log2标准化: log2(sparse.gbm.C + 1 )
对行 求 和: Matrix::rowSums(sparse.gbm)
对行求均值: Matrix::rowMeans(sparse.gbm)
对列 求 和: Matrix::colSums(sparse.gbm)
对行求均值: Matrix::colMeans(sparse.gbm)
5、读写出稀疏矩阵
- 读取
readMM( "matrix.mtx") - 写出
writeMM(sparse.gbm.C, "matrix.mtx");system("gzip --keep matrix.mtx")
参考资料
Working with a sparse matrix in R - Kamil Slowikowski (slowkow.com)
Instructions for using Python SciPy sparse matrix | Develop Paper
R convert matrix or data frame to sparseMatrix - Stack Overflow
网友评论