刘小泽写于19.7.18
之前接触过scRNA的Seurat包 2.X版本，后来2019.4.16正式升级到了3，虽然有一些函数进行了调整和拆分，但总体思路上还是变化不大的，这次就来探索一下。因为这个是个大包，所以需要写几篇才能系统学完
前言
还是先说下什么是Seurat吧，它是一个分析单细胞转录组数据的R包，从QC到最后的分群、基因鉴定，可以探索单细胞的异质性，也可以整合不同的单细胞数据类型。
进入主页(https://satijalab.org/seurat/)，就看到开心的官宣：我们终于升级到3.0啦！作者酝酿了差不多一年，终于推出了正式版，直接使用install.packages('Seurat')安装即可。
简单说几点在版本2上的更新：
-
增加扩展了单细胞整合的方法。虽说世界上没有两片相同的树叶，但也能找到不同事物之间的关联。它就是通过找不同细胞类型之间的
anchors(可以理解为一个锚，将两个不同细胞类型固定在一起)，构建"和谐"的细胞群（可参考：Stimulated vs. Control PBMCs），或者在不同实验之间传递信息(可参考：Multiple Dataset Integration and Label Transfer) -
改进标准化方法。提到一种新方法：
sctransform，它有一个最显著的特点就是：Compared to standard log-normalization, sctransform effectively removes technically-driven variation while preserving biological heterogeneity.
并且强调拒绝了：We do not apply heuristic steps, such as log-transformation, pseudocount addition, or z-scoring; do not assume a constant size or global ‘scaling factor’ for single cells
它的教程在：https://satijalab.org/seurat/v3.0/sctransform_vignette.html
-
重构了Seurat对象，更新点在于多数据整合，并且可以轻松切换multi-modal data
Easily switch between RNA, protein, cell hashing, batch-corrected / integrated, or imputed data.
当然，一个好的开发团队绝不会对之前的作品置之不理，Seurat同时保留了Seurat2的使用，包括：
- 它的安装（https://satijalab.org/seurat/install.html#previous）；
- 使用说明(在官网右下角帮助按钮，选择版本2即可；这一点和10X做的很像，兼顾了多种版本)；
- 另外版本3提供了
UpgradeSeuratObject函数，可以分析之前的版本2数据，使用方法很简单：updated_seurat_object = UpdateSeuratObject(object = old_seurat_object)就会返回一个版本3的数据结构； - 提供了大量的版本2和版本3的对比
cheatsheet(https://satijalab.org/seurat/essential_commands.html)
标准的Seurat分析流程
虽然发展到了3.0阶段，但核心与2.0也类似，基本都是通过跑一个流程去探索数据的分群信息。一般包括：数据normalization、找有差异的feature(应该也不限于基因，因为版本3将以前版本2的很多带gene的函数都改成了feature ，可能是为了体现新版本的包容性更强)、数据scale、PCA、构建SNN（shared-nearest-neighbors）图、利用一个叫modularity optimizer的算法分群，最后用t-SNE可视化
大体上操作就是：
pbmc.counts <- Read10X(data.dir = "~/Downloads/pbmc3k/filtered_gene_bc_matrices/hg19/")
pbmc <- CreateSeuratObject(counts = pbmc.counts)
pbmc <- NormalizeData(object = pbmc)
pbmc <- FindVariableFeatures(object = pbmc)
pbmc <- ScaleData(object = pbmc)
pbmc <- RunPCA(object = pbmc)
pbmc <- FindNeighbors(object = pbmc)
pbmc <- FindClusters(object = pbmc)
pbmc <- RunTSNE(object = pbmc)
DimPlot(object = pbmc, reduction = "tsne")
关于normalization、scale、standardization：https://kharshit.github.io/blog/2018/03/23/scaling-vs-normalization
normalization：调整数据符合正态分布（Normal distribution/Gaussian distribution/ bell curve），因为一些下游统计或者机器学习算法需要用到正态分布的数据，比如 t-tests, ANOVAs, linear regression, linear discriminant analysis (LDA) and Gaussian Naive Bayes等
scale：调整数据范围，例如落在[0，1]中
standardization：(also called z-score normalization) ：将数据变成均值为0，标准差为1的分布，在 SVM, logistics regression and neural networks中应用较多
Seurat3.0 新的交互式操作
相比2.0，重新设计的函数确实人性化不少，例如以前取列名(也就是细胞样本名)需要用object@cell.names ，也就是对象的操作；而新版貌似将操作"降了一级"，只要colnames(x = object)就好，很像数据框或矩阵的操作
此外还有：
# 得到细胞名(两种方法)
colnames(x = pbmc)
Cells(object = pbmc)
> identical(Cells(PBMC),colnames(PBMC))
[1] TRUE
# 得到基因名
rownames(x = pbmc)
# 得到数量信息
ncol(x = pbmc)
nrow(x = pbmc)
# Get cell identity classes
Idents(object = pbmc)
levels(x = pbmc)
# Stash cell identity classes
pbmc[["old.ident"]] <- Idents(object = pbmc)
pbmc <- StashIdent(object = pbmc, save.name = "old.ident")
# Set identity classes
Idents(object = pbmc) <- "CD4 T cells"
Idents(object = pbmc, cells = 1:10) <- "CD4 T cells"
# Set identity classes to an existing column in meta data
Idents(object = pbmc, cells = 1:10) <- "orig.ident"
Idents(object = pbmc) <- "orig.ident"
# Rename identity classes
pbmc <- RenameIdents(object = pbmc, `CD4 T cells` = "T Helper cells")
## 取子集
# 可以在identity class的基础上取
subset(x = pbmc, idents = "B cells")
subset(x = pbmc, idents = c("CD4 T cells", "CD8 T cells"), invert = TRUE)
# 可以根据基因/feature的表达量取
subset(x = pbmc, subset = MS4A1 > 3)
# 可以设定取值范围
subset(x = pbmc, subset = MS4A1 > 3 & PC1 > 5)
subset(x = pbmc, subset = MS4A1 > 3, idents = "B cells")
# 根据对象的 meta data取
subset(x = pbmc, subset = orig.ident == "Replicate1")
# 每个identity class取一定数量的细胞
subset(x = pbmc, downsample = 100)
## 组合
# 两个对象
merge(x = pbmc1, y = pbmc2)
# 组合多于两个seurat对象
merge(x = pbmc1, y = list(pbmc2, pbmc3))
发现新版的设计让对象的操作变得更灵活，更贴近我们日常对数据框的操作了
访问数据
访问meta信息
# 查看存在object@meta.data中的meta信息
pbmc[[]]
# 提取特定的meta信息
pbmc$nCount_RNA
pbmc[[c("percent.mito", "nFeature_RNA")]]
# 添加meta信息(版本2中使用AddMetaData，当然在版本3也能用这种方式)，但这里又是类似于数据框的操作，易于理解
random_group_labels <- sample(x = c("g1", "g2"), size = ncol(x = pbmc), replace = TRUE)
pbmc$groups <- random_group_labels
访问表达量信息
# Retrieve or set data in an expression matrix ('counts', 'data', and 'scale.data')
GetAssayData(object = pbmc, slot = "counts")
pbmc <- SetAssayData(object = pbmc, slot = "scale.data", new.data = new.data)
访问降维后的细胞/feature信息
Embeddings(object = PBMC, reduction = "pca") #细胞
Loadings(object = PBMC, reduction = "pca") #feature
Loadings(object = pbmc, reduction = "pca", projected = TRUE)
任意访问
# 同时选取多种数据
# FetchData can pull anything from expression matrices, cell embeddings, or metadata
FetchData(object = pbmc, vars = c("PC_1", "percent.mito", "MS4A1"))
Seurat的教程是相当丰富的，不用担心学不会，一步步跟着走，跑一遍流程，心里就有底了。
接下来会有几篇推送，跟着官网https://satijalab.org/seurat/vignettes.html教程走一遍，看看这个包到底能给我们带来什么惊喜
欢迎关注我们的公众号～_～
我们是两个农转生信的小硕，打造生信星球，想让它成为一个不拽术语、通俗易懂的生信知识平台。需要帮助或提出意见请后台留言或发送邮件到jieandze1314@gmail.com
Welcome to our bioinfoplanet!
网友评论