R语言-稀疏矩阵对象格式学习-重点理解稀疏矩阵对象的重构

稀疏矩阵几乎产生于所有的大型科学工程计算领域，记录样本特征值的稠密矩阵中很多记录值都是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( )方法转换 稠密矩阵 $\Rightarrow$ 稀疏矩阵

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 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( )方法转换 稀疏矩阵 $\Rightarrow$ 稠密矩阵

 >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