這一篇文章是回答知識星球中一位星友的提問，她的電腦內(nèi)存有限，無法直接使用所有數(shù)據(jù)，只能分析部分數(shù)據(jù)。

數(shù)據(jù)來源: https://content.cruk.cam.ac.uk/jmlab/atlas_data.tar.gz 解壓縮之后，得到下面數(shù)據(jù)

數(shù)據(jù)清單

其中raw_counts.mtx是以稀疏矩陣格式存放的表達量數(shù)據(jù)，文件為6.5G, 用普通的文本編輯器無法打開，我們可以用Linux命令行的less查看數(shù)據(jù)存放形式

稀疏矩陣

顯然這種格式并不是給人類閱讀的，它存放的是非零數(shù)據(jù)的位置及其具體數(shù)值。當然，我們也不需要讀懂，只需要R語言或者其他編程語言能夠加載即可。

R語言的Matrix包的readMM函數(shù)就能夠讀取該文件

mt <- Matrix::readMM("raw_counts.mtx")
dim(mt)
# [1]  29452 139331
# 行為基因，列為細胞

這一步時間非常的久，我差不多花了10分鐘時間。同時占用內(nèi)存也非?？捎^，直接占用了8G左右的內(nèi)存，不到16G內(nèi)存的電腦可能根本無法讀取。

format(object.size(mt), units = "Mb")
# "7377.8 Mb"

稀疏矩陣其實和普通矩陣看起來差不多，除了在顯示的時候用.來表示0.

查看部分信息

還有一點就是，對于這種量級的數(shù)據(jù)，我們無法使用R自帶的as.data.frame或者as.matrix將其轉(zhuǎn)成普通的數(shù)據(jù)庫或者矩陣，它會直接報錯。因此我也不建議對其進行數(shù)據(jù)轉(zhuǎn)換。

我們發(fā)現(xiàn)這里的矩陣并沒有行名和列名，這部分信息需要額外從其他文件中讀取

bc <- read.table("barcodes.tsv")
genes <- read.table("genes.tsv", sep = "\t")
dim(bc)
#[1] 139331      1
dim(genes)
#[1] 29452     2

不難發(fā)現(xiàn)barcode的行數(shù)等于矩陣的列數(shù)， gene的行數(shù)等于矩陣的行數(shù)， 也就是說矩陣的列是細胞，行是基因。

row.names(mt) <- genes$V1
colnames(mt) <- bc$V1

查看部分信息

建議：將此處得到matrix保存為Rds格式，方便后續(xù)加載

saveRDS(mt, "raw_matrix.Rds")

接下來就是根據(jù)元信息來提取對應(yīng)的細胞，我們以提取"Mesenchyme"細胞為例進行講解

meta.info <- read.table("meta.tab",
                        sep = "\t", header = TRUE)

cell.info <- meta.info[meta.info$celltype == "Mesenchyme", "cell"]
cell.info <- cell.info[!is.na(cell.info)]

mt.sml <- mt[, cell.info]
format(object.size(mt.sml), units = "Mb")
# "280.9 Mb"

代碼的核心邏輯為提取出對應(yīng)行的細胞名，然后根據(jù)細胞名提取矩陣中的對應(yīng)列。

過濾后的細胞就可以用作后續(xù)分析。不過在開始分析之前，讓我們先把原始的矩陣給刪掉，因為它實在是太占用內(nèi)存了。

rm(mt); gc()

除了用元信息進行過濾外，你還可以通過隨機抽樣，從原始數(shù)據(jù)中抽出部分細胞，這樣子也能夠在內(nèi)存吃緊的情況進行后續(xù)分析。