在之前推文python scipy 稀疏矩陣詳解中，詳細(xì)介紹了常用稀疏矩陣原理及在python中的使用經(jīng)驗。本篇推文聚焦R語言稀疏矩陣格式及其在單細(xì)胞多組學(xué)(scRNA, scATAC)中的應(yīng)用。

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，是一種按列壓縮的稀疏矩陣格式，由三個一維數(shù)組indptr, indices, data組成。dgCMatrxi是R包Matrix中標(biāo)準(zhǔn)的格式，也是最常見的格式，單細(xì)胞多組學(xué)稀疏數(shù)據(jù)存儲一般采用該格式。

dgRMatrix

即Sparse Compressed, Row-oriented Numeric Matrices，類比于Python中的csr_matrix，和dgCMatrxi相反，該格式稀疏矩陣按行壓縮，原理上同樣由三個一維數(shù)組indptr, indices, data組成。dgRMatrxi稀疏矩陣在單細(xì)胞組學(xué)分析中一般不顯示使用。所以接下來的介紹以dgCMatrxi為主，dgRMatrxi的實現(xiàn)方法實際是類似的。

R稀疏矩陣常用操作方法

稀疏矩陣構(gòu)建

######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內(nèi)置函數(shù)，用于組建常規(guī)矩陣。 mat3 <- as(mat3, "dgCMatrix") # dgCMatrix #3 x 3 sparse Matrix of class "dgCMatrix" #[1,] . 2 . #[2,] . . . #[3,] 1 3 .方法四：new("dgCMatrix", ...)，...處填寫dgCMatrix類的屬性具體值 #雖然R是統(tǒng)計語言，但R也可以面向?qū)ο缶幊獭4對象是標(biāo)準(zhǔn)的面向?qū)ο髮崿F(xiàn)方式，有專門的類定義函數(shù)setClass()和類的實例化函數(shù)new()。"dgCMatrix" 本質(zhì)上就是一個類，所以我們可以直接用new函數(shù)來實例化這個類創(chuàng)建對象。 indptr <- c(2L, 0L, 2L) #必為整數(shù) indices <- c(0L, 1L, 3L, 3L) #必為整數(shù) data <- c(1, 2, 3) mat4 <- new('dgCMatrix', i=indptr, p=indices, x=data, Dim=c(3L,3L)) #Dim參數(shù)必為整數(shù)######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") #和構(gòu)建dgC矩陣是一樣的，不贅述方法四：new("dgTMatrix", ...) #和構(gòu)建dgC矩陣是一樣的，不贅述

稀疏矩陣格式轉(zhuǎn)換

#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") # 創(chuàng)建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包通常用來操作稀疏矩陣，而內(nèi)置函數(shù)matrix則用來操作常規(guī)矩陣。 #[1] FALSE

稀疏矩陣讀寫

#通常稀疏矩陣的讀取和保存是通過Matrix包實現(xiàn)的。 Matrix::readMM #讀取稀疏矩陣 Matrix::writeMM #保存稀疏矩陣

稀疏矩陣單細(xì)胞多組學(xué)中的應(yīng)用

Cellranger

以下是Cellranger生成的結(jié)果，其中matrix.mtx文件是dgTMatrix格式，所以如果下游分析軟件要求其他格式的稀疏矩陣則需要有轉(zhuǎn)換的步驟。

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

Seurat

使用Seurat和Scanpy分析單細(xì)胞數(shù)據(jù)一個很大的不同在于：Annadata對象中每一行為一個細(xì)胞，而Seurat對象則相反。Seurat中單細(xì)胞稀疏數(shù)據(jù)存儲采用dgCMatrix；而Cellranger輸出到文件的稀疏存儲方式是dgTMatrix格式，所以用Seurat分析Cellranger輸出的數(shù)據(jù)必然要先做稀疏矩陣格式的轉(zhuǎn)換，而 Seurat::Read10X函數(shù)的核心實現(xiàn)就是這個， Seurat::Read10X函數(shù)會生成帶有行列名的dgCMatrix。當(dāng)然你也可以不用這個函數(shù)，即自己構(gòu)建稀疏矩陣然后利用CreateSeuratObject函數(shù)生成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') #將稀疏矩陣轉(zhuǎn)換為Seurat默認(rèn)的dgC格式 colnames(mtx) = row.names(barcodes) rownames(mtx) = row.names(genes) pbmc = CreateSeuratObject(counts = mtx, assay = 'RNA', meta.data = barcodes)

Signac

對于 Cellranger 數(shù)據(jù)的導(dǎo)入可以使用 Seurat::Read10X 和 Seurat::Read10X_h5函數(shù)。但當(dāng)我們從公共數(shù)據(jù)庫中下載數(shù)據(jù)時，很多情況下我們只能下載到cellranger-atac 產(chǎn)生的三個文件：barcodes.tsv、matrix.mtx和peaks.bed。那么這種情況下，我們要么強(qiáng)行嘗試 Seurat::Read10X函數(shù)，要么明晰原理自己實現(xiàn)一個讀取方式。

如果你直接用 Seurat::Read10X讀取這三個文件會發(fā)現(xiàn)代碼報錯，原因有兩個：

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

Seurat::Read10X 函數(shù)需要讀入不重復(fù)的barcodes和features(gene name or peak ID)，但是peaks.bed并沒有提供unique的peak ID。

因此，方法一，修改peaks.bed文件內(nèi)容，然后使用Seurat::Read10X函數(shù)讀取

# peaks.bed文件格式，沒有unique的peak ID， 需要對其進(jìn)行轉(zhuǎn)換。 chr1 9757 10652 chr1 86735 87645 chr1 180601 181169 #向peaks.bed文件中添加peak ID，構(gòu)建的方法同上；然后再修改文件名；最后用Read10X函數(shù)讀取。 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='-') #構(gòu)建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函數(shù)讀取。 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")

方法二，自己讀取數(shù)據(jù)，構(gòu)建對象

#Read10X_h5函數(shù)本質(zhì)上就是生成有行名和列名的dgCMatrix，所以如果我們沒有h5文件，我們可以直接自己構(gòu)建這個矩陣。 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='-') #構(gòu)建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生成的稀疏矩陣， 用此對象繼續(xù)下游分析。 chrom_assay <- CreateChromatinAssay(counts = mtx, sep = c(":", "-"), genome = 'hg19',# fragments = '../vignette_data/atac_v1_pbmc_10k_fragments.tsv.gz', #沒有這個文件可以不提供，這樣后面不能做fragments相關(guān)的分析min.cells = 10,min.features = 200) pbmc <- CreateSeuratObject(counts = chrom_assay, assay = "peaks")

參考資料

R Matrix

Signac

公眾號-編程自修室：稀疏矩陣詳細(xì)解析與單細(xì)胞數(shù)據(jù)讀取

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。