Seurat Weekly NO.4 || 高效数据管理

2020-12-03 本文已影响0人周运来就是我

在这里，和国际同行一起学习单细胞数据分析。

单细胞转录组数据分析框架是成熟的，所谓降维聚类必知必会嘛。这一点也不影响数据分析过程的非线性过程，我们知道，在数据科学中降维聚类均属于探索性数据分析手段。探索时贪婪，验证时谨慎。在处理单细胞转录组数据的时候，往往时螺旋式的，不时会有反复。这个过程就需要我们有数据管理的策略，好在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为我们提供了toolsSlot。也就是不管我们引入的函数也好，自己写的函数也好，均可与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秒"

牺牲存储换取时间，数十倍的差别，数据稍微大点，这就够喝杯咖啡的了。

参考：
https://github.com/satijalab/seurat/wiki
https://www.cnblogs.com/maoerbao/p/11694958.html

Seurat Weekly NO.4 || 高效数据管理

猜你喜欢

热点阅读