【原】Seurat V5|当单细胞进入百万细胞时代，BPCell 给出一种“解”决参考

生信补给站 2023-12-26 发布于北京

展开全文

现在单细胞的样本量和数据量都在逐渐增多，大有百万之势，如果用电脑分析就会显的吃力。这时候Seurat V5通过调用BPCells 也许能给出一个解决方案。

一 R包，数据准备

1，准备BPCells等R包

如果因网络问题出现报错的话，可以在github上下载BPCells的zip文件，然后本地安装。之后缺少什么包就安装什么包。

#devtools::install_github(c("bnprks/BPCells"))library(Seurat)library(BPCells)library(ggplot2)# needs to be set for large dataset analysisoptions(future.globals.maxSize = 1e9)

2，准备数据

使用https://cf./samples/cell-exp/1.3.0/1M_neurons/1M_neurons_filtered_gene_bc_matrices_h5.h5中下载的h5数据进行示例，细胞数130万。

##数据读取brain.data <- open_matrix_10x_hdf5(path = "1M_neurons_filtered_gene_bc_matrices_h5.h5")##输出表达矩阵到本地write_matrix_dir(mat = brain.data,dir = "brain_counts")

会在当前目录下根据write_matrix_dir生产一个对应名字的文件夹（brain_counts），内含有很多文件用于后续读取。

#查看数据brain.mat <- open_matrix_dir(dir = "brain_counts")brain.mat

直接读取10X的三个文件会报错。

二标准分析

1，创建Seurat 对象

根据前面读取的矩阵构建object对象，进行HVG之前的标准分析

# Create Seurat Objectpbmc <- CreateSeuratObject(counts = mat_raw)pbmc
##标准化数据brain <- NormalizeData(brain)##识别高变基因brain <- FindVariableFeatures(brain)

head(brain)orig.ident nCount_RNA nFeature_RNAAAACCTGAGATAGGAG-1 SeuratProject       4046         1807AAACCTGAGCGGCTTC-1 SeuratProject       2087         1249AAACCTGAGGAATCGC-1 SeuratProject       4654         2206AAACCTGAGGACACCA-1 SeuratProject       3193         1655AAACCTGAGGCCCGTT-1 SeuratProject       8444         3326AAACCTGAGTCCGGTC-1 SeuratProject      11178         3866

2，抽取细胞

使用SketchData函数抽取细胞，ncells参数指定抽取多少细胞到内存中进行分析，注意指定assay名字，方便后续分析。其他细胞仍然可以硬盘访问。

##抓取2w细胞载入到内存中进行分析brain <- SketchData(object = brain,ncells = 20000,method = "LeverageScore",sketched.assay = "sketch")brain
An object of class Seurat 55996 features across 1306127 samples within 2 assays Active assay: sketch (27998 features, 2000 variable features)2 layers present: counts, data1 other assay present: RNA
dim(brain@assays$sketch$counts)#[1] 27998 20000

可以看到已完成2W细胞的抽取，接下来进行后续分析。

3，Sketch细胞标准分析

将默认对象切换到sketch（内存中的2w细胞），然后进行标准分析

# 默认对象切换到sketch（内存中的2w细胞）DefaultAssay(brain) <- "sketch"brain <- FindVariableFeatures(brain)brain <- ScaleData(brain)brain <- RunPCA(brain)brain <- FindNeighbors(brain, dims = 1:50)brain <- FindClusters(brain, resolution = 2)brain <- RunUMAP(brain, dims = 1:50, return.model = T)brain

除了上述sketch抽取的2W外，其他的cluster信息（sketch_snn_res.2）均为NA 。

（2）基于2W细胞进行可视化

##绘制umap图（2w细胞）p1 <- DimPlot(brain, label = T, label.size = 3, reduction = "umap") + NoLegend()p1p2 <- FeaturePlot(object = brain,                  features = c("ENSMUSG00000036256",           "ENSMUSG00000037984",           "ENSMUSG00000023391",           "ENSMUSG00000026787"),                  ncol = 2)p2

4 扩展至全部细胞

使用ProjectData函数将sketch得到的cluster和reduction信息扩展至全部百万细胞。

扩展后的降维，聚类结果存放在含有full的结果中，如pca.full， full.umap和cluster_full 。

##扩展（从2w细胞扩展到130w细胞）brain <- ProjectData(object = brain,assay = "RNA",full.reduction = "pca.full",sketched.assay = "sketch",sketched.reduction = "pca",umap.model = "umap",dims = 1:50,refdata = list(cluster_full = "seurat_clusters"))

新增的cluster_full即为根据sketch的cluster预测出来的全部细胞的cluster结果，cluster_full.score为预测的得分。

更多细节可以通过View(TransferSketchLabels)查看

（2）基于全部细胞可视化

##默认对象切换到RNA（全部数据集）DefaultAssay(brain) <- "RNA"p3 <- DimPlot(brain, label = T,               label.size = 3,               reduction = "full.umap", group.by = "cluster_full",               alpha = 0.1) + NoLegend()p1 + p3

（3）以两个基因为例，比较sketch 和全部的结果

DefaultAssay(brain) <- "sketch"x1 <- FeaturePlot(brain, features = c("ENSMUSG00000036256", "ENSMUSG00000036887") )DefaultAssay(brain) <- "RNA"x2 <- FeaturePlot(brain, features = c("ENSMUSG00000036256", "ENSMUSG00000036887") )p4 <- x1 / x2p4