- Seurat Weekly NO.0 || 开刊词
- Seurat Weekly NO.1 || 到底分多少个群是合适的？！
- Seurat Weekly NO.2 || 我该如何取子集
- Seurat Weekly NO.3 || 直接用Seurat画fig2
在这里，和国际同行一起学习单细胞数据分析。
单细胞转录组数据分析框架是成熟的，所谓降维聚类必知必会嘛。这一点也不影响数据分析过程的非线性过程，我们知道，在数据科学中降维聚类均属于探索性数据分析手段。探索时贪婪，验证时谨慎。在处理单细胞转录组数据的时候，往往时螺旋式的，不时会有反复。这个过程就需要我们有数据管理的策略，好在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
网友评论