R语言-稀疏矩阵对象格式介绍&重构方法

背景：稀疏矩阵几乎产生于所有的大型科学工程计算领域，记录样本特征值的稠密矩阵中很多记录值都是0，由于0不携带信息，因此耗费空间存储0元素是很浪费资源的行为。而且很多计算只对非零元素进行操作，将特征矩阵构建成稀疏矩阵，可以很容易的索引到非零元素，所以基于稀疏矩阵的数据运算，可是实现更低的资源占用和更快的计算速度。

在单细胞领域，稀疏矩阵对于处理 scRNA-seq 表达谱数据是非常必要的，构建分析对象的时候它们消耗更低的内存。本文重点介绍 R语言 稀疏矩阵对象格式， 稀疏矩阵与稠密矩阵的相互转换。

1、稀疏矩阵的格式类别

• 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## 按行向量方向排列的非零元素 值 （从上到下）

2、R 语言操作稀疏矩阵的软件包-Matrix

Matrix 是R语言中用于处理矩阵和稀疏矩阵的一个重要包，可以提高数据处理的效率和速度，特别是在处理大型数据集时。。它提供了一种新的矩阵类，称为“稀疏矩阵”，可以更高效地处理大型数据集。Matrix包中的函数可以用于创建、操作和计算矩阵和稀疏矩阵，包括矩阵的转置、乘法、加法、求逆等操作。Matrix包还提供了一些特殊的矩阵类型，如对称矩阵、三角矩阵和对角矩阵等，这些矩阵类型在某些情况下可以提高计算效率。Matrix包还支持多线程计算，可以利用多核CPU提高计算速度。参考： Matrix package - RDocumentation

install.packages("Matrix") ; library(Matrix)

2.1 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
  

2.2 重构sparseMatrix

2.2.1 分解sparseMatrix再 重构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 #矩阵的维度
    )

$\color{red}{注意：}$由于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

2.2.2 从long-data数据框 重构sparseMatrix

参考： R - how to convert long-data dataframe to sparse matrix - Stack Overflow

### PrePare Matrix
pacman::p_load(Matrix,dplyr)
mt_wider <- data.frame(matrix(rpois(12,1),nrow = 3))
colnames(mt_wider) <- paste0("c",1:4); rownames(mt_wider) <- paste0("r",1:3)

### Format wider DT to longer DT
mt_long <- mt_wider %>% tibble::rownames_to_column("col_name") %>% reshape2::melt(id = c("col_name")) %>% 
    filter(value != 0) %>% rename("row_name"=2) %>% 
    mutate(across(1:2,.fns = as.factor)) %>% tibble()
####################################################################################
> mt_long
# A tibble: 9 × 3
  col_name row_name value
  <fct>    <fct>    <int>
1 r1       c1           2
2 r2       c1           1
3 r3       c1           3
4 r1       c2           1
5 r2       c2           1
6 r1       c3           3
7 r2       c3           2
8 r3       c3           1
9 r1       c4           2
####################################################################################

### Format longer DT to SparseMTX
mt_sparse <- with(mt_long, sparseMatrix(i=as.numeric(col_name), 
                                        j=as.numeric(row_name), 
                                        x=as.numeric(value),
                                        dimnames=list(levels(col_name), 
                                                      levels(row_name))))

### Check for Restoration
> as.matrix(mt_sparse) == mt_wider
     c1   c2   c3   c4
r1 TRUE TRUE TRUE TRUE
r2 TRUE TRUE TRUE TRUE
r3 TRUE TRUE TRUE TRUE

2.3 稀疏矩阵的一些运算

log2标准化: log2(sparse.gbm.C + 1 )
对行 求 和: Matrix::rowSums(sparse.gbm)
对行求均值: Matrix::rowMeans(sparse.gbm)
对列 求 和: Matrix::colSums(sparse.gbm)
对行求均值: Matrix::colMeans(sparse.gbm)

2.4 读写出稀疏矩阵

• 读取
  readMM( "matrix.mtx")
• 写出
  writeMM(sparse.gbm.C, "matrix.mtx");system("gzip --keep matrix.mtx")

Reference

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
R - how to convert long-data dataframe to sparse matrix - Stack Overflow