面对几十万细胞,首先要区分“保留完整表达数据”与“每一步都在完整数据上计算”。Seurat的磁盘矩阵与sketch流程提供一种可选路径:完整counts留在磁盘,选定细胞用于部分探索,再把学到的表示或标签投影到全量数据。它不是把大任务自动变成完全无内存需求的操作。
下面演示当前公开工作流,要求Seurat、SeuratObject和BPCells兼容,保存辅助函数还需要fs等依赖。例子以已经确认是RNA基因表达计数的10x HDF5为输入;含ADT、其他模态或不同特征命名的文件应先按模态核对,不能直接把所有行当作RNA基因。
使用新的工作目录,保留原HDF5不动。BPCells负责磁盘矩阵的读写,Seurat保存对象结构和路径。创建之后核对assay类、特征和细胞名,不调用as.matrix,也不根据一个很小的RDS大小判断完整数据只占很少空间。
library(Seurat)
library(SeuratObject)
library(BPCells)
dir.create("analysis_bundle", showWarnings=FALSE)
bp_dir <- "analysis_bundle/counts_bp"
if (dir.exists(bp_dir)) stop("请使用新的输出目录,避免覆盖已有矩阵")
mat <- open_matrix_10x_hdf5("filtered_feature_bc_matrix.h5")
write_matrix_dir(mat=mat, dir=bp_dir)
counts <- open_matrix_dir(bp_dir)
obj <- CreateSeuratObject(counts=counts)
stopifnot(inherits(obj[["RNA"]], "Assay5"),
identical(Layers(obj[["RNA"]]), "counts"))
print(list(dim=dim(obj[["RNA"]]), class=class(obj[["RNA"]]$counts)))
writeLines(capture.output(sessionInfo()), "analysis_bundle/session.txt")
下面的细胞数、维度和分辨率仅为演示参数,要用目标数据检验。示例单counts层意味着ncells对应本层抽取量;多层时当前接口允许按层指定数量,单个整数也可能按每层解释,不能误以为总数永远等于该整数。
DefaultAssay(obj) <- "RNA"
obj <- NormalizeData(obj)
obj <- FindVariableFeatures(obj)
obj <- SketchData(obj, assay="RNA", ncells=min(10000L, ncol(obj)),
method="LeverageScore", sketched.assay="sketch", seed=2026L)
DefaultAssay(obj) <- "sketch"
obj <- FindVariableFeatures(obj)
obj <- ScaleData(obj, features=VariableFeatures(obj))
npcs <- min(30L, length(VariableFeatures(obj))-1L, ncol(obj[["sketch"]])-1L)
stopifnot(npcs>=2L, ncol(obj[["sketch"]])>50L)
obj <- RunPCA(obj, npcs=npcs)
obj <- FindNeighbors(obj, dims=seq_len(npcs))
obj <- FindClusters(obj, resolution=0.5)
obj <- RunUMAP(obj, dims=seq_len(npcs), return.model=TRUE)| 对象或步骤 | 主要作用 | 必须核对的边界 |
|---|---|---|
| RNA磁盘counts | 保留全部原始计数访问 | 后端目录是数据的一部分 |
| sketch assay | 支持一部分探索和建模 | 不是完整数据的同义替身 |
| PCA与UMAP模型 | 在抽取细胞上学习表示 | 参数与稀有群体表现需检查 |
| 全量投影 | 扩展所选表示或标签 | 不等于所有步骤全量重算 |
| 保存与迁移 | 保留对象与后端依赖 | 不能只复制一个RDS文件 |
用与sketch一致的维度和UMAP模型投影,并分别命名全量结果。得到全量标签后,检查细胞覆盖和已知标记,特别关注稀有群体是否被过度合并。sketch有自己的抽样目标,标签投影也有近似性,不能凭图形好看保证每个稀有群体都被正确发现。
obj <- ProjectData(obj, assay="RNA", sketched.assay="sketch",
sketched.reduction="pca", full.reduction="pca.full",
umap.model="umap", dims=seq_len(npcs),
normalization.method="LogNormalize",
refdata=list(cluster_full="seurat_clusters"))
DefaultAssay(obj) <- "RNA"
stopifnot(nrow(Embeddings(obj, "pca.full"))==ncol(obj[["RNA"]]))
print(table(is.na(obj$cluster_full)))
SaveSeuratRds(obj, file="analysis_bundle/analysis.rds",
move=FALSE, relative=TRUE)
check <- LoadSeuratRds("analysis_bundle/analysis.rds")
print(dim(check[["RNA"]]))
此处后端原本建在bundle内,保存时显式不移动已有数据,并使用相对路径。将完整目录复制到新位置后,还应在新R会话重新载入、抽查少量值和必要降维。不要在验证前清理原始输入或后端目录;对象能被读取但外部路径失效,后续访问仍会失败。
HDF5文件能够被打开,只说明读取路径可用。仍需按实验情况确认低质量细胞、空液滴及特征命名等问题,并保存质量筛选依据。不能为了让对象更小就随意删掉难处理的样本,也不能把没有做质量检查的示例流程直接当作最终科学结果。
在全量投影前,先在sketch中查看与研究问题相关的标记和样本覆盖。如果某个关键群体根本没有合适的参考细胞,投影不会自动提供新的证据。对需要深入分析的群体,应返回原始细胞集合建立适用的后续方案,并说明该方案的计算范围。
大规模绘图也可能产生额外时间与内存需求。先用适当范围核查表达和标签,再决定是否输出全量图;绘图慢、后端路径失效和模型计算不足,应分别排查。最终交接时,保留目录清单、版本与关键检查结果,才能让新的服务器验证同一套对象,而不是只收到一个可以打开外壳的RDS。目标账号还应能读取全部后端目录,权限验收也属于交接检查的一部分。
分别测量写盘、归一化、sketch、投影和保存的内存与I/O峰值。后续算法如果不支持该矩阵类或会物化数据,仍可能需要较大内存。确认必要步骤和资源缺口后,再评估适合的大内存CPU、工作盘与租期,服务器租赁的环境部署和数据保留范围需逐项确认。