稀疏矩陣幾乎產(chǎn)生于所有的大型科學(xué)工程計算領(lǐng)域，記錄樣本特征值的稠密矩陣中很多記錄值都是0，由于0不攜帶信息，因此耗費空間存儲0元素是很浪費資源的行為。而且很多計算只對非零元素進行操作，將特征矩陣構(gòu)建成稀疏矩陣，可以很容易的索引到非零元素，所以基于稀疏矩陣的數(shù)據(jù)運算，可是實現(xiàn)更低的資源占用和更快的計算速度。

在單細胞領(lǐng)域，稀疏矩陣對于處理 scRNA-seq 表達譜數(shù)據(jù)是非常必要的，構(gòu)建分析對象的時候它們消耗更低的內(nèi)存。本文重點介紹 R語言 稀疏矩陣對象格式，稀疏矩陣的重構(gòu)， 稀疏矩陣與稠密矩陣的相互轉(zhuǎn)換。

三種類型的稀疏矩陣對象結(jié)構(gòu)

首先，稀疏矩陣的一般格式主要有三種：

• dgTMatrix ::: 創(chuàng)建未壓縮的 三元組 稀疏矩陣 - 坐標格式的稀疏矩陣 - coo
  coo格式的稀疏矩陣以三元組(i,j,x)的形式記錄非零條目
  • i ## 非零元素的行號【0-based索引值】
  • j ## 非零元素的列號【0-based索引值】
  • x ## 按列向量方向排列的非零元素 值（從左到右）
• dgCMatrix ::: 創(chuàng)建列壓縮的 三元組 稀疏矩陣 - 壓縮系數(shù)矩陣 - csc - 更快的列檢索速度
  csc格式的稀疏矩陣以三元組(i,p,x)的形式記錄非零條目
  • i ## 非零元素的行號【0-based索引值】
  • p ## 從0開始, 記錄值是對每列非零元素的個數(shù)的累加（從左到右）
  • x ## 按列向量方向排列的非零元素 值 （從左到右）
• dgRMatrix ::: 創(chuàng)建行壓縮的 三元組 稀疏矩陣 - csr - 更快的行檢索速度
  • p ## 從0開始, 記錄值是對每行非零元素的個數(shù)的累加（從上到下）
  • j ## 非零元素的列號【0-based索引值】
  • x ## 按行向量方向排列的非零元素 值 （從上到下）

1、R 語言 Matrix包

????? Matrix 包為擴展基本矩陣數(shù)據(jù)類型的密集和稀疏矩陣提供了一組 S4 類和應(yīng)用于這些類對象的各種函數(shù)和運算符方法。>Matrix package - RDocumentation

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

2、R 語言 稠密矩陣 與 稀疏矩陣 的相互轉(zhuǎn)換

構(gòu)建一個測試用稠密矩陣：

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("矩陣非零元素個數(shù) ", sum( colSums(counts != 0)) ," 個") )
[1] "矩陣非零元素個數(shù) 12 個"
> is(counts, 'sparseMatrix') #判斷為非 稀疏矩陣 對象
[1] FALSE

• 通過 as( ) 方法轉(zhuǎn)換 稠密矩陣 稀疏矩陣

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

創(chuàng)建列壓縮的稀疏矩陣(csc_matrix)還可以用方法sparse.gbm <- Matrix::Matrix(counts, sparse = T)，但是該方法在大型矩陣中較為低效。
三種稀疏矩陣對象在R語言中的結(jié)構(gòu)str(sparse.gbm.*)

數(shù)據(jù)結(jié)構(gòu)

關(guān)于壓縮格式的稀疏矩陣對象中數(shù)據(jù)結(jié)構(gòu)的 sparse.gbm.C@p 向量，該向量不是非零元素的行/列索引編號，它的值是對 數(shù)據(jù)壓縮方向上 每個向量內(nèi)非零元素個數(shù)的累加(如創(chuàng)建列壓縮稀疏矩陣時，壓縮方向為按列 向左到右 移動)。p向量的記錄方式相比i/j記錄非零元素索引號，有助于進一步減少稀疏矩陣對象存儲大小。

p的計算

#通過 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( ) 方法轉(zhuǎn)換 稀疏矩陣 稠密矩陣

 >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、重構(gòu)sparseMatrix

sparse.gbm.C.new <- Matrix::sparseMatrix(
    i = sparse.gbm.C@i + 1 ,  # +1 轉(zhuǎn)換 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中數(shù)據(jù)對象（如data.frame,vector）的索引值都是 【1-based】 的，所以在通過Matrix::sparseMatrix重構(gòu)稀疏矩陣的時候，i/j的索引值要加1轉(zhuǎn)換為【1-based】

除了按S4的方法從對象中取三元組信息，通過Matrix::summary()方法可以快速獲取稀疏矩陣S4對象的（i,j,x)三元組，并用于重構(gòu) 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