在之前推文python scipy 稀疏矩阵详解中，详细介绍了常用稀疏矩阵原理及在python中的使用经验。本篇推文聚焦R语言稀疏矩阵格式及其在单细胞多组学(scRNA, scATAC)中的应用。

R稀疏矩阵

dgTMatrix

即Sparse matrices in triplet form，该格式类比于python中的coo_matrix，是最简单的稀疏矩阵存储方式，采用三元组(row, col, data)(或称为ijv format)的形式来存储矩阵中非零元素的信息。

dgCMatrix

即Compressed, sparse, column-oriented numeric matrices，类比于Python中的csc_matrix，是一种按列压缩的稀疏矩阵格式，由三个一维数组indptr, indices, data组成。dgCMatrxi是R包Matrix中标准的格式，也是最常见的格式，单细胞多组学稀疏数据存储一般采用该格式。

dgRMatrix

即Sparse Compressed, Row-oriented Numeric Matrices，类比于Python中的csr_matrix，和dgCMatrxi相反，该格式稀疏矩阵按行压缩，原理上同样由三个一维数组indptr, indices, data组成。dgRMatrxi稀疏矩阵在单细胞组学分析中一般不显示使用。所以接下来的介绍以dgCMatrxi为主，dgRMatrxi的实现方法实际是类似的。

R稀疏矩阵常用操作方法

稀疏矩阵构建

######dgCMatrix######
方法一：Matrix(*, sparse = TRUE)
#Matrix包
mat <- Matrix(data = 0L, nrow=3201, ncol = 818, sparse = TRUE) #dgCMatrix
print(object.size(mat), unit="KB")
#4.7 Kb
dim(mat)
#[1] 3201  818

方法二：sparseMatrix()
#Matrix包
i <- c(1,3:8)
j <- c(2,9,6:10)
x <- 7 * (1:7)
mat2 <- sparseMatrix(i, j, x = x) # dgCMatrix

方法三：as(*, "dgCMatrix")
mat3 <- matrix(data=c(0,0,1,2,0,3,0,0,0), nrow=3, ncol=3) #matrix内置函数，用于组建常规矩阵。
mat3 <- as(mat3, "dgCMatrix") # dgCMatrix
#3 x 3 sparse Matrix of class "dgCMatrix"      
#[1,] . 2 .
#[2,] . . .
#[3,] 1 3 .

方法四：new("dgCMatrix", ...)，...处填写dgCMatrix类的属性具体值
#虽然R是统计语言，但R也可以面向对象编程。S4对象是标准的面向对象实现方式，有专门的类定义函数setClass()和类的实例化函数new()。"dgCMatrix" 本质上就是一个类，所以我们可以直接用new函数来实例化这个类创建对象。
indptr <-  c(2L, 0L, 2L) #必为整数
indices <- c(0L, 1L, 3L, 3L) #必为整数
data <- c(1, 2, 3)
mat4 <- new('dgCMatrix', i=indptr, p=indices, x=data, Dim=c(3L,3L))  #Dim参数必为整数



######dgTMatrix######
方法一：sparseMatrix()
#Matrix包
i <- c(1,3:8)
j <- c(2,9,6:10)
x <- 7 * (1:7)
mat1 <- sparseMatrix(i, j, x = x, repr = "T") 

方法二：spMatrix()
#Matrix包
mat2 <- spMatrix(10,20, i = c(1,3:8),
                     j = c(2,9,6:10),
                     x = 7 * (1:7))

方法三：as(*, "dgTMatrix") #和构建dgC矩阵是一样的，不赘述

方法四：new("dgTMatrix", ...) #和构建dgC矩阵是一样的，不赘述

稀疏矩阵格式转换

#dense matrix to sparse matrix
sparse_matrix <- as(dense_matrix, "dgCMatrix")
sparse_matrix <- as(dense_matrix, "dgTMatrix")
sparse_matrix <- as(dense_matrix, 'sparseMatrix')  # dgCMatrix
sparse_matrix <- as(dense_matrix, 'Matrix')  # dgCMatrix

#sparse matrix to dense matrix
dense_matrix <- as.matrix(sparse_matrix)
dense_matrix <- as(sparse_matrix, 'matrix')

#dgCMatrix to dgTMatrix
dgC_matrix <- as(dgT_matrix, 'dgCMatrix')
dgT_matrix <- as(dgC_matrix, 'dgTMatrix')

其实还有一个R包SparseM包能非常方便的操作稀疏矩阵，这个R包中稀疏矩阵的名称和python中是一样的，即csc、csr和coo稀疏矩阵。感兴趣的读者可自行了解。

稀疏矩阵判断

有时候需要判断矩阵是否为稀疏矩阵以及是哪种稀疏矩阵，可用如下命令：

i <- c(1,3:8)
j <- c(2,9,6:10)
x <- 7 * (1:7)
mat <- sparseMatrix(i, j, x = x, repr = "T") # 创建dgTMatrix
is(mat, 'sparseMatrix') 
#[1] TRUE
is(mat, 'dgTMatrix') 
#[1] TRUE
is(mat, 'dgCMatrix')
#[1] FALSE
is(mat, 'Matrix') 
#[1] TRUE
is(mat, 'matrix') #需要注意的是Matrix包通常用来操作稀疏矩阵，而内置函数matrix则用来操作常规矩阵。
#[1] FALSE

稀疏矩阵读写

#通常稀疏矩阵的读取和保存是通过Matrix包实现的。
Matrix::readMM #读取稀疏矩阵
Matrix::writeMM #保存稀疏矩阵

稀疏矩阵单细胞多组学中的应用

Cellranger

以下是Cellranger生成的结果，其中matrix.mtx文件是dgTMatrix格式，所以如果下游分析软件要求其他格式的稀疏矩阵则需要有转换的步骤。

$ tree filtered_feature_bc_matrix/
filtered_feature_bc_matrix/
├── barcodes.tsv.gz
├── features.tsv.gz  
└── matrix.mtx.gz  dgT格式

Seurat

使用Seurat和Scanpy分析单细胞数据一个很大的不同在于：Annadata对象中每一行为一个细胞，而Seurat对象则相反。Seurat中单细胞稀疏数据存储采用dgCMatrix；而Cellranger输出到文件的稀疏存储方式是dgTMatrix格式，所以用Seurat分析Cellranger输出的数据必然要先做稀疏矩阵格式的转换，而 Seurat::Read10X函数的核心实现就是这个， Seurat::Read10X函数会生成带有行列名的dgCMatrix。当然你也可以不用这个函数，即自己构建稀疏矩阵然后利用CreateSeuratObject函数生成Seurat对象。参考代码如下：

## construct RNA Seurat object
dir.10x = './filtered_feature_bc_matrix'
genes = read.delim(paste0(dir.10x, 'genes.tsv'), row.names = 1)
barcodes = read.delim(paste0(dir.10x, 'barcodes.tsv'), row.names = 1)
mtx = Matrix::readMM(paste0(dir.10x, 'matrix.mtx'))
mtx = as.matrix(mtx)
mtx = as(mtx, 'dgCMatrix') #将稀疏矩阵转换为Seurat默认的dgC格式
colnames(mtx) = row.names(barcodes)
rownames(mtx) = row.names(genes)
pbmc = CreateSeuratObject(counts = mtx, assay = 'RNA', meta.data = barcodes)

Signac

对于 Cellranger 数据的导入可以使用 Seurat::Read10X 和 Seurat::Read10X_h5函数。但当我们从公共数据库中下载数据时，很多情况下我们只能下载到cellranger-atac 产生的三个文件：barcodes.tsv、matrix.mtx和peaks.bed。那么这种情况下，我们要么强行尝试 Seurat::Read10X函数，要么明晰原理自己实现一个读取方式。

如果你直接用 Seurat::Read10X读取这三个文件会发现代码报错，原因有两个：

1. 文件名不对。因为 Seurat::Read10X 默认读入的三个文件名是barcodes.tsv.gz、features.tsv.gz和matrix.mtx.gz，而cellranger-atac输出的文件名是有差异的。
2. Seurat::Read10X 函数需要读入不重复的barcodes和features(gene name or peak ID)，但是peaks.bed并没有提供unique的peak ID。

因此，方法一，修改peaks.bed文件内容，然后使用Seurat::Read10X函数读取

# peaks.bed文件格式，没有unique的peak ID， 需要对其进行转换。
chr1	9757	10652
chr1	86735	87645
chr1	180601	181169

#向peaks.bed文件中添加peak ID，构建的方法同上；然后再修改文件名；最后用Read10X函数读取。
data.dir <- './filtered_peak_bc_matrix'
features <- read.csv('./filtered_peak_bc_matrix/peaks.bed',  header = FALSE, row.names = NULL, sep = '\t')
row.names(features) <- paste(paste(features$V1, features$V2, sep=':'), features$V3, sep='-')   #构建peak ID
features <- features['V1']  #单独取出peak ID这一列保存
write.table(features, './filtered_peak_bc_matrix/features.tsv', sep = '\t', row.names=F, col.names=F, quote = F)
#接下来再修改其他的文件名，即可用Read10X函数读取。
counts <- Read10X(data.dir,gene.column = 1, cell.column = 1)
chrom_assay <- CreateChromatinAssay(counts = counts, 
                                    sep = c(":", "-"), genome = 'hg19',
  																	# fragments = '../vignette_data/atac_v1_pbmc_10k_fragments.tsv.gz', #没有这个文件可以不提供，这样后面不能画track图
  																	min.cells = 10,
  																	min.features = 200)
pbmc <- CreateSeuratObject(counts = chrom_assay, assay = "peaks")

方法二，自己读取数据，构建对象

#Read10X_h5函数本质上就是生成有行名和列名的dgCMatrix，所以如果我们没有h5文件，我们可以直接自己构建这个矩阵。
dir.10x = './filtered_feature_bc_matrix'
features = read.csv(paste0(dir.10x, 'peaks.bed'), header = FALSE, row.names = NULL, sep = '\t') 
row.names(features) <- paste(paste(features$V1, features$V2, sep=':'), features$V3, sep='-')  #构建peak ID
barcodes = read.delim(paste0(dir.10x, 'barcodes.tsv'), row.names = 1)
mtx = Matrix::readMM(paste0(dir.10x, 'matrix.mtx'))
colnames(mtx) = row.names(barcodes)
mtx = as(mtx, 'dgCMatrix') #此处生成的稀疏矩阵即等价于Read10X_h5生成的稀疏矩阵， 用此对象继续下游分析。
chrom_assay <- CreateChromatinAssay(counts = mtx, 
                                    sep = c(":", "-"), genome = 'hg19',
  																	# fragments = '../vignette_data/atac_v1_pbmc_10k_fragments.tsv.gz', #没有这个文件可以不提供，这样后面不能做fragments相关的分析
  																	min.cells = 10,
  																	min.features = 200)
pbmc <- CreateSeuratObject(counts = chrom_assay, assay = "peaks")

参考资料

R Matrix

Signac

公众号-编程自修室：稀疏矩阵详细解析与单细胞数据读取