「单细胞转录组系列」如何从稀疏矩阵中随机提取数据进行分析

「单细胞转录组系列」如何从稀疏矩阵中提取部分数据进行分析写了根据细胞的元信息来提取数据，然后有人问我，如何进行随机抽样，这篇文章就来介绍下如何随机提取细胞。

这次的数据集在https://data.humancellatlas.org/explore/projects/cc95ff89-2e68-4a08-a234-480eca21ce79/expression-matrices 大约有100K细胞。我下载的其中mtx格式，下载链接如下

https://data.humancellatlas.org/project-assets/project-matrices/cc95ff89-2e68-4a08-a234-480eca21ce79.homo_sapiens.mtx.zip

经过解压缩之后，我们会得到5个文件，其中matrix.mtx.gz存放表达量矩阵，压缩后有2.2G

3.5M Nov  7  2019 barcodes.tsv.gz # barcode信息
109M Nov  7  2019 cells.tsv.gz    # 细胞元信息   
1.1M Nov  7  2019 features.tsv.gz # 基因的元信息
1.1M Nov  7  2019 genes.tsv.gz    # 基因的元信息, 同features, 有header
2.2G Nov  7  2019 matrix.mtx.gz   # 表达量信息

我们用R读取其中cells.tsv.gz和matrix.mtx.gz, 加载表达量矩阵和细胞信息

mtx <- readMM("e7448a34-b33d-41de-b422-4c09bfeba96b.mtx/matrix.mtx.gz")
dim(mtx)
# 58347 782859
paste0("Memory Size = ", round(object.size(mtx) / 10^6, 3), " MB")
# "Memory Size = 9908.946 MB"
cells <- data.table::fread("e7448a34-b33d-41de-b422-4c09bfeba96b.mtx/cells.tsv.gz")
dim(cells)
# 782859     44

一共有78K细胞，占用了接近10G的内存。为了降低电脑分析压力，我们可以提取部分的细胞进行分析。

第一种抽样方法，就是用sample无放回的抽取

cell_id <- sample(x = seq(1, nrow(cells)), 10000)
mtx_sml <- mtx[, cell_id]
paste0("Memory Size = ", round(object.size(mtx_sml) / 10^6, 3), " MB")
#  "Memory Size = 125.365 MB
cell_sml <- cells[cell_id, ]

第二种抽样方法，我们可以根据细胞的元信息，分层抽样，例如donor_organism.human_specific.ethnicity.ontology_label记录着细胞着来源，于是我们就可以根据不同人群进行抽样。

# 获取细胞编号
cell_ids <- row.names(cells)
# 获取细胞来源
cell_source <- factor(cells$donor_organism.human_specific.ethnicity.ontology_label)
# 根据细胞来源换分细胞，然后抽样
cell_ids_list <- lapply(split(cell_ids, cell_source), function(x){ sample(x, 3000) })
# 提取细胞编号
cell_id <- unlist(cell_ids_list)
# 提取细胞
cell_sml <- cells[cell_id, ]
mtx_sml <- mtx[, cell_id]

我们可以统计抽样后不同类型的细胞占比，发现细胞占比是一样的。

table(cell_sml$donor_organism.human_specific.ethnicity.ontology_label)
#                 African American            Asian         European 
#            3000             3000             3000             3000

当然，如果原来细胞元信息中有细胞类型注释信息，还可以按照细胞分群抽样。