Seurat Weekly NO.4 || 高效数据管理
天子呼来不上船,
自称臣是菜鸟团。
在这里,和国际同行一起学习单细胞数据分析。
单细胞转录组数据分析框架是成熟的,所谓降维聚类必知必会
嘛。这一点也不影响数据分析的非线性过程,我们知道,在数据科学中降维聚类均属于探索性数据分析手段。探索时贪婪,验证时谨慎。在处理单细胞转录组数据的时候,往往时螺旋式的,不时会有反复。这个就需要我们有数据管理的策略,好在Seurat的S4结构为我们的数据管理做好了准备。
记录分析过程
数据分析的时候,会需要看不同参数下的数据表现,比如分群的resolution
,差异分析的不同方法(或不同分组)。这时候我们如果能够记录下来每次分析的命令/时间等,在复盘或重现分析的时候就方便很多。Seurat为我们提供了commands
Slot。如我们可以这样查看之前的分析命令:
library(Seurat)
library(SeuratData)
pbmc3k.final -> pbmc
pbmc@commands$RunUMAP.RNA.pca
Command: RunUMAP(pbmc3k.final, dims = 1:10)
Time: 2020-04-30 12:55:01
dims : 1 2 3 4 5 6 7 8 9 10
reduction : pca
assay : RNA
umap.method : uwot
n.neighbors : 30
n.components : 2
metric : cosine
learning.rate : 1
min.dist : 0.3
spread : 1
set.op.mix.ratio : 1
local.connectivity : 1
repulsion.strength : 1
negative.sample.rate : 5
uwot.sgd : FALSE
seed.use : 42
angular.rp.forest : FALSE
verbose : TRUE
reduction.name : umap
reduction.key : UMAP_
Seurat的主要分析函数都内置了LogSeuratCommand
函数,在完成计算之后,把计算参数记录在pbmc@commands
中。我们当然也可以自己往这里面写入提示自己的信息:
pbmc@commands$菜鸟团数据分析 <- stringr::str_glue(paste0(Sys.time(),"\ninitiation: \n你的名字"))
pbmc@commands$菜鸟团数据分析
2020-12-03 20:37:28
initiation:
你的名字
自己写的函数
进阶阶段肯定是在自己的单细胞转录分析函数中加入LogSeuratCommand
了,这样我们对数据做的每一次分析参数都会保留下来了。那么,问题来了,我们的函数可不可以写在Seurat内部呢?
显然是可以的,Seurat为我们提供了tools
Slot。也就是不管我们引入的函数也好,自己写的函数也好,均可与Seurat合为一体,当然,函数名要注意一下,一如Tools函数的帮助文档里备注的一样:
For developers: set tool data using Tool<-. Tool<- will automatically set the name of the tool to the function that called Tool<-, so each function gets one entry in the tools list and cannot overwrite another function’s entry. The automatic naming will also remove any method identifiers (eg. RunPCA.Seurat will become RunPCA); please plan accordingly.
存储数据
我们知道Seurat为单细胞转录组分析过程存下了counts/data/scale.data三张表达谱,还还有一系列的降维的坐标等。如果我们有临床信息还可以存在meta.data里,但是要是还有其他信息呢?比如差异基因列表。这时候,可以存在misc
Slot里面,省的每次重新计算或者省去写出读入的操作。
Idents(pbmc) <- "seurat_clusters"
pbmc@misc$clustersmk <- FindAllMarkers(pbmc)
Idents(pbmc) <- "seurat_annotations"
pbmc@misc$celltypesmk <- FindAllMarkers(pbmc)
如有感兴趣的基因集也可以存在其中,经典的marker存在这里不是很方便吗?
pbmc@misc$chemokines <- grep("^CXC|CCL|CCR|CX3|XCL|XCR", rownames(pbmc),value = T)
快速读取
随着单细胞的通量越来越大,我们的Seurat对象也越来越大(另一方面是Seurat存的数据越来越多)。每天来到办公室第一件事就是读入一个Seurat对象,如何才能快速读取呢?介绍saveRDS的一个参数。
?saveRDS
compress
a logical specifying whether saving to a named file is to use “gzip” compression, or one of “gzip”, “bzip2” or “xz” to indicate the type of compression to be used. Ignored if file is a conne
t1=proc.time()
saveRDS(pbmc,'pbmc1.rds',compress = F)
t2=proc.time()
t=t2-t1
print(paste0('执行时间:',t[3][[1]],'秒'))
[1] "执行时间:1.83000000000015秒"
t1=proc.time()
saveRDS(pbmc,'pbmc.rds')
t2=proc.time()
t=t2-t1
print(paste0('执行时间:',t[3][[1]],'秒'))
[1] "执行时间:27.6600000000001秒"
牺牲存储换取时间,数十倍的差别,数据稍微大点,这就够喝杯咖啡的了。
以上,是我们基于Seurat结构给出的数据管理的解决方案。数据分析过程的记录同实验记录是一样重要的,从第一次读入数据到文章online都尽量做到:可重复/可交付。
文末友情推荐
参考:
https://github.com/satijalab/seurat/wiki
https://www.cnblogs.com/maoerbao/p/11694958.html