參考：https://satijalab.org/seurat/articles/essential_commands.html

????

1. Seurat 標(biāo)準(zhǔn)工作流程

讀入數(shù)據(jù) -> 構(gòu)建Seurat對(duì)象 -> 歸一化 -> 尋找特征基因 -> 標(biāo)準(zhǔn)化 -> PCA降維 -> 計(jì)算KNN最近鄰 -> 分群 -> 非線性降維UMAP和TSNE -> 可視化降維結(jié)果

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")

2. Seurat對(duì)象交互

獲取細(xì)胞和基因名，并計(jì)數(shù)

colnames(x = pbmc)
rownames(x = pbmc)
ncol(x = pbmc)
nrow(x = pbmc)

設(shè)置identity classes（active.ident）

# 查看目前所用的identity（也就是查看目前用于細(xì)胞分類的metadata）并計(jì)數(shù)
Idents(object = pbmc)
levels(x = pbmc)

# 存放細(xì)胞identity classes
pbmc[["old.ident"]] <- Idents(object = pbmc) #將"old.ident"存放在pbmc的metadata中

# 設(shè)置identity classes
Idents(object = pbmc) <- "CD4 T cells"  #將"CD4 T cells"設(shè)置為pbmc的 "active.ident"
Idents(object = pbmc, cells = 1:10) <- "CD4 T cells" #將前10個(gè)細(xì)胞的ident設(shè)置為"CD4 T cells"

# 將meta data中已有的分類設(shè)置為identity classes
Idents(object = pbmc, cells = 1:10) <- "orig.ident" #將前10個(gè)細(xì)胞的ident設(shè)置為"orig.ident"
Idents(object = pbmc) <- "orig.ident" #將"orig.ident"設(shè)置為pbmc的 "active.ident"

# 重命名 identity classes
pbmc <- RenameIdents(object = pbmc, `CD4 T cells` = "T Helper cells")

從seurat對(duì)象中取子集??

# 根據(jù)identity class取子集（also see ?SubsetData）
subset(x = pbmc, idents = "B cells")
subset(x = pbmc, idents = c("CD4 T cells", "CD8 T cells"), invert = TRUE)

# 根據(jù)某個(gè)基因/特征的表達(dá)水平取子集
subset(x = pbmc, subset = MS4A1 > 3)

# 根據(jù)組合的標(biāo)準(zhǔn)取子集
subset(x = pbmc, subset = MS4A1 > 3 & PC1 > 5)
subset(x = pbmc, subset = MS4A1 > 3, idents = "B cells")

# 根據(jù)某個(gè)meta data中的某一分組來(lái)取子集
subset(x = pbmc, subset = orig.ident == "Replicate1")

# Downsample the number of cells per identity class
subset(x = pbmc, downsample = 100)

多個(gè)Seurat對(duì)象的合并

# Merge兩個(gè)Seurat對(duì)象
merge(x = pbmc1, y = pbmc2)
# Merge多個(gè)Seurat對(duì)象
merge(x = pbmc1, y = list(pbmc2, pbmc3))

3. 數(shù)據(jù)訪問(wèn)

#查看metadata數(shù)據(jù)庫(kù)框, metadata存儲(chǔ)在object@meta.data
pbmc[[]]

# 檢索metadata中的特定指標(biāo)
pbmc$nCount_RNA
pbmc[[c("percent.mito", "nFeature_RNA")]]

# 添加metadata, （見?AddMetaData）
random_group_labels <- sample(x = c("g1", "g2"), size = ncol(x = pbmc), replace = TRUE)
pbmc$groups <- random_group_labels

# 檢索或設(shè)置表達(dá)矩陣數(shù)據(jù)，包括原始表達(dá)矩陣、標(biāo)準(zhǔn)化的矩陣和降維后的矩陣 ('counts', 'data', and 'scale.data')
GetAssayData(object = pbmc, slot = "counts")[1:5,1:5]
count.data <- GetAssayData(object = pbmc[["RNA"]], slot = "counts")
count.data <- as.matrix(x = count.data + 1)
new.seurat.object <- SetAssayData(object = pbmc, slot = "counts", new.data = count.data, assay = "RNA"）
#?? 使用GetAssayData函數(shù)可以從Seurat對(duì)象訪問(wèn)數(shù)據(jù)。
#?? 可以使用SetAssayData將數(shù)據(jù)添加到counts，data或scale.data插槽中。新數(shù)據(jù)必須具有與當(dāng)前數(shù)據(jù)相同順序的相同細(xì)胞。添加到counts'或data`中的數(shù)據(jù)必須具有與當(dāng)前數(shù)據(jù)相同的features。

# Get cell embeddings and feature loadings
Embeddings(object = pbmc, reduction = "pca") #檢索每個(gè)細(xì)胞的PC矩陣
Loadings(object = pbmc, reduction = "pca") #檢索每個(gè)基因的PC矩陣
Loadings(object = pbmc, reduction = "pca", projected = TRUE)

# ??FetchData函數(shù)可以從expression matrices, cell embeddings或metadata中取出任何值
head(FetchData(object = pbmc, vars = c("PC_1", "percent.mt", "MS4A1")))
#                        PC_1 percent.mt    MS4A1
# AAACATACAACCAC-1 -4.7296855  3.0177759 0.000000
# AAACATTGAGCTAC-1 -0.5174029  3.7935958 2.583047
# AAACATTGATCAGC-1 -3.1891063  0.8897363 0.000000
# AAACCGTGCTTCCG-1 12.7933021  1.7430845 0.000000
# AAACCGTGTATGCG-1 -3.1288078  1.2244898 0.000000
# AAACGCACTGGTAC-1 -3.1088963  1.6643551 0.000000

4. 可視化

Seurat的繪圖基于ggplot2
詳見Seurat繪圖函數(shù)總結(jié)

# 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圖中的預(yù)制個(gè)性化主題

# 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)

5. Multi-Assay Features

With Seurat, you can easily switch between different assays at the single cell level (such as ADT counts from CITE-seq, or integrated/batch-corrected data). Most functions now take an assay parameter, but you can set a Default Assay to avoid repetitive statements.

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")

6. 其他好用的函數(shù)

HVFInfo函數(shù)從Assay對(duì)象中提取特征均值和離散度。可變特征向量可以通過(guò)VariableFeatures函數(shù)提取。VariableFeatures也可以設(shè)置可變特征向量。

# HVFInfo pulls mean, dispersion, and dispersion scaled
# Useful for viewing the results of FindVariableFeatures
> head(x = HVFInfo(object = pbmc))
#                      mean    variance variance.standardized
# AL627309.1    0.003411676 0.003401325             0.9330441
# AP006222.2    0.001137225 0.001136363             0.9924937
# RP11-206L10.2 0.001895375 0.001892500             0.9627290
# RP11-206L10.9 0.001137225 0.001136363             0.9924937
# LINC00115     0.006823351 0.006779363             0.9062135
# NOC2L         0.107278241 0.159514698             0.7849309

# VariableFeatures both accesses and sets the vector of variable features
> head(x = VariableFeatures(object = pbmc))
# [1] "PPBP"   "LYZ"    "S100A9" "IGLL5"  "GNLY"   "FTL" 
# Set variable features example

可以通過(guò)Stdev找到Seurat對(duì)象中存儲(chǔ)的DimReduc的標(biāo)準(zhǔn)差向量。

Stdev(object = pbmc, reduction.use = 'pca')
# Warning: The following arguments are not used: reduction.use
#  [1] 7.098420 4.495493 3.872592 3.748859 3.171755 2.545292 2.068137 1.945133
#  [9] 1.847375 1.834689 1.820439 1.788429 1.779215 1.757395 1.751558 1.746384
# [17] 1.732116 1.729550 1.722981 1.720541 1.717668 1.715182 1.710343 1.705735
# [25] 1.705204 1.702498 1.700713 1.695816 1.694965 1.690976 1.688388 1.681423
# [33] 1.679596 1.678656 1.674677 1.674487 1.670084 1.667662 1.666598 1.664131
# [41] 1.660317 1.657907 1.656963 1.654456 1.654188 1.649514 1.647186 1.645571
# [49] 1.642952 1.641145

可以通過(guò)以下方法找到Seurat類：

library(Seurat)
utils::methods(class = 'Seurat')