Seurat4.0系列教程6：常用命令
Seurat 标准流程
标准 Seurat 工作流采用原始的单细胞表达数据，旨在数据中查找clusters。此过程包括数据标准化和高变基因选择、数据归一化、高变基因的PCA、共享近邻图形的构建以及使用模块优化进行聚类。最后，我们使用 t-SNE 在二维空间中可视化我们的clusters。
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")
Seurat对象交互
自 Seurat v4.0 以来，对 Seurat 对象进行了改进，并增加了用户交互的新方法。还为常见任务引入简单的功能，例如取子集和合并。
# Get cell and feature names, and total numbers
colnames(x = pbmc)
Cells(object = pbmc)
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")# Subset Seurat object based on identity class, also see ?SubsetData
subset(x = pbmc, idents = "B cells")
subset(x = pbmc, idents = c("CD4 T cells", "CD8 T cells"), invert = TRUE)
# Subset on the expression level of a gene/feature
subset(x = pbmc, subset = MS4A1 > 3)
# Subset on a combination of criteria
subset(x = pbmc, subset = MS4A1 > 3 & PC1 > 5)
subset(x = pbmc, subset = MS4A1 > 3, idents = "B cells")
# Subset on a value in the object meta data
subset(x = pbmc, subset = orig.ident == "Replicate1")
# Downsample the number of cells per identity class
subset(x = pbmc, downsample = 100)# Merge two Seurat objects
merge(x = pbmc1, y = pbmc2)
# Merge more than two Seurat objects
merge(x = pbmc1, y = list(pbmc2, pbmc3))
数据访问
在 Seurat 中访问数据很简单，使用明确定义的取子集代码可以快速查找所需的数据。
# View metadata data frame, stored in object@meta.data
pbmc[[]]
# Retrieve specific values from the metadata
pbmc$nCount_RNA
pbmc[[c("percent.mito", "nFeature_RNA")]]
# Add metadata, see ?AddMetaData
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)# Get cell embeddings and feature loadings
Embeddings(object = pbmc, reduction = "pca")
Loadings(object = pbmc, reduction = "pca")
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的可视化
默认情况下，所有绘图功能都将返回 ggplot2 绘图，从而允许使用 ggplot2 轻松定制。
# Dimensional reduction plot for PCA or tSNE
DimPlot(object = pbmc, reduction = "tsne")
DimPlot(object = pbmc, reduction = "pca")
# Dimensional reduction plot, with cells colored by a quantitative feature
FeaturePlot(object = pbmc, features = "MS4A1")
# Scatter plot across single cells, replaces GenePlot
FeatureScatter(object = pbmc, feature1 = "MS4A1", feature2 = "PC_1")
FeatureScatter(object = pbmc, feature1 = "MS4A1", feature2 = "CD3D")
# Scatter plot across individual features, repleaces CellPlot
CellScatter(object = pbmc, cell1 = "AGTCTACTAGGGTG", cell2 = "CACAGATGGTTTCT")
VariableFeaturePlot(object = pbmc)
# Violin and Ridge plots
VlnPlot(object = pbmc, features = c("LYZ", "CCL5", "IL32"))
RidgePlot(object = pbmc, feature = c("LYZ", "CCL5", "IL32"))
# Heatmaps
DoHeatmap(object = pbmc, features = heatmap_markers)
DimHeatmap(object = pbmc, reduction = "pca", cells = 200)
# New things to try! Note that plotting functions now return ggplot2 objects, so you can add themes, titles, and options
# onto them
VlnPlot(object = pbmc, features = "MS4A1", split.by = "groups")
DotPlot(object = pbmc, features = c("LYZ", "CCL5", "IL32"), split.by = "groups")
FeaturePlot(object = pbmc, features = c("MS4A1", "CD79A"), blend = TRUE)
DimPlot(object = pbmc) + DarkTheme()
DimPlot(object = pbmc) + labs(title = "2,700 PBMCs clustered using Seurat and viewed\non a two-dimensional tSNE")
Seurat 提供了许多预制的主题，可添加到 ggplot2 中，以便快速定制
| 主题 | 功能 |
|---|---|
DarkTheme |
设置带有白色文本的黑色背景 |
FontSize |
为绘图的各个元素设置字体大小 |
NoAxes |
删除轴和轴文本 |
NoLegend |
删除所有图例元素 |
RestoreLegend |
删除后恢复图例 |
RotatedAxis |
旋转 x 轴标签 |
# Plotting helper functions work with ggplot2-based scatter plots, such as DimPlot, FeaturePlot, CellScatter, and
# FeatureScatter
plot <- DimPlot(object = pbmc) + NoLegend()
# HoverLocator replaces the former `do.hover` argument It can also show extra data throught the `information` argument,
# designed to work smoothly with FetchData
HoverLocator(plot = plot, information = FetchData(object = pbmc, vars = c("ident", "PC_1", "nFeature_RNA")))
# FeatureLocator replaces the former `do.identify`
select.cells <- FeatureLocator(plot = plot)
# Label points on a ggplot object
LabelPoints(plot = plot, points = TopCells(object = pbmc[["pca"]]), repel = TRUE)
不同assay之间切换
使用 Seurat，您可以轻松地在单细胞水平的不同assays 之间切换（例如来自 CITE-seq 的 ADT 计数，或整合/批次校正的数据）。大多数功能现在需要assays 参数，但可以设置默认assay以避免重复语句。
cbmc <- CreateSeuratObject(counts = cbmc.rna)
# Add ADT data
cbmc[["ADT"]] <- CreateAssayObject(counts = cbmc.adt)
# Run analyses by specifying the assay to use
NormalizeData(object = cbmc, assay = "RNA")
NormalizeData(object = cbmc, assay = "ADT", method = "CLR")
# Retrieve and set the default assay
DefaultAssay(object = cbmc)
DefaultAssay(object = cbmc) <- "ADT"
DefaultAssay(object = cbmc)
# Pull feature expression from both assays by using keys
FetchData(object = cbmc, vars = c("rna_CD3E", "adt_CD3"))
# Plot data from multiple assays using keys
FeatureScatter(object = cbmc, feature1 = "rna_CD3E", feature2 = "adt_CD3")
Seurat v2.X与v4.X命令的区别
| Seurat v2.X | Seurat v4.X |
|---|---|
| object@data | GetAssayData(object = object)[1] |
| object@raw.data | GetAssayData(object = object, slot = "counts") |
| object@scale.data | GetAssayData(object = object, slot = "scale.data") |
| object@cell.names | colnames(x = object)[2] |
| rownames(x = object@data)[3] | rownames(x = object)[4] |
| object@var.genes | VariableFeatures(object = object)[5] |
| object@hvg.info | HVFInfo(object = object)[6] |
| object@assays$assay.name | object[["assay.name"]] |
| object@dr$pca | object[["pca"]] |
| GetCellEmbeddings(object = object, reduction.type = "pca") | Embeddings(object = object, reduction = "pca") |
| GetGeneLoadings(object = object, reduction.type = "pca") | Loadings(object = object, reduction = "pca") |
| AddMetaData(object = object, metadata = vector, col.name = "name") | object$name <- vector |
| object@meta.data$name | object$name |
| object@idents | Idents(object = object)[7] |
| SetIdent(object = object, ident.use = "new.idents") | Idents(object = object) <- "new.idents" |
| SetIdent(object = object, cells.use = 1:10, ident.use = "new.idents") | Idents(object = object, cells = 1:10) <- "new.idents" |
| StashIdent(object = object, save.name = "saved.idents") | object$saved.idents <- Idents(object = object) |
| levels(x = object@idents)[8] | levels(x = object)[9] |
| RenameIdent(object = object, old.ident.name = "old.ident", new.ident.name = "new.ident") | RenameIdents(object = object, "old.ident" = "new.ident") |
| WhichCells(object = object, ident = "ident.keep") | WhichCells(object = object, idents = "ident.keep") |
| WhichCells(object = object, ident.remove = "ident.remove") | WhichCells(object = object, idents = "ident.remove", invert = TRUE) |
| WhichCells(object = object, max.cells.per.ident = 500)[10] | WhichCells(object = object, downsample = 500)[11] |
| WhichCells(object = object, subset.name = "name", low.threshold = low, high.threshold = high) | WhichCells(object = object, expression = name > low & name < high)[12] |
| FilterCells(object = object, subset.names = "name", low.threshold = low, high.threshold = high) | subset(x = object, subset = name > low & name < high)[13] |
| SubsetData(object = object, subset.name = "name", low.threshold = low, high.threshold = high) | subset(x = object, subset = name > low & name < high)[14] |
| MergeSeurat(object1 = object1, object2 = object2) | merge(x = object1, y = object2)[15] |
文中链接
GetAssayData(object = object): https://satijalab.org/seurat/reference/AssayData.html
[2]
colnames(x = object): https://rdrr.io/r/base/colnames.html
[3]
rownames(x = object@data): https://rdrr.io/r/base/colnames.html
[4]
rownames(x = object): https://rdrr.io/r/base/colnames.html
[5]
VariableFeatures(object = object): https://satijalab.org/seurat/reference/VariableFeatures.html
[6]
HVFInfo(object = object): https://satijalab.org/seurat/reference/VariableFeatures.html
[7]
Idents(object = object): https://satijalab.org/seurat/reference/Idents.html
[8]
levels(x = object@idents): https://rdrr.io/r/base/levels.html
[9]
levels(x = object): https://rdrr.io/r/base/levels.html
[10]
WhichCells(object = object, max.cells.per.ident = 500): https://satijalab.org/seurat/reference/WhichCells.html
[11]
WhichCells(object = object, downsample = 500): https://satijalab.org/seurat/reference/WhichCells.html
[12]
WhichCells(object = object, expression = name > low & name < high): https://satijalab.org/seurat/reference/WhichCells.html
[13]
subset(x = object, subset = name > low & name < high): https://rdrr.io/r/base/subset.html
[14]
subset(x = object, subset = name > low & name < high): https://rdrr.io/r/base/subset.html
[15]
merge(x = object1, y = object2): https://rdrr.io/r/base/merge.html