scRNA-单细胞Seurat包升级，换汤不换药

Jul 18, 2019 5 min read cnposts

刘小泽写于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.

    说明在： Multimodal vignette

当然，一个好的开发团队绝不会对之前的作品置之不理，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")
# 但是好像访问速度比原来的object@assays$RNA要慢一点
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教程走一遍，看看这个包到底能给我们带来什么惊喜

​

scRNA
Yunze Liu
Yunze Liu
Bioinformatics Sharer

Co-founder of Bioinfoplanet(生信星球)

Next
Previous

Related