在之前推文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讀取這三個文件會發現代碼報錯，原因有兩個：

文件名不對。因為 Seurat::Read10X 默認讀入的三個文件名是barcodes.tsv.gz、features.tsv.gz和matrix.mtx.gz，而cellranger-atac輸出的文件名是有差異的。

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

公眾號-編程自修室：稀疏矩陣詳細解析與單細胞數據讀取

總結

以上是生活随笔為你收集整理的R语言稀疏矩阵详解的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。