Seurat4.0系列教程1：标准流程
分享是一种态度
时代的洪流奔涌而至，单细胞技术也从旧时王谢堂前燕，飞入寻常百姓家。雪崩的时候，没有一片雪花是无辜的，你我也从素不相识，到被一起卷入单细胞天地。R语言和Seurat已以势如破竹之势进入4.0时代，天问一号探测器已进入火星乌托邦平原了，你还不会单细胞吗？那么为了不被时代抛弃的太远，跟着我们一起通过学习seurat系列教程入门单细胞吧。首先我们学习经典的标准流程，这个教程小编当初入门时候曾经花1000购买过，也曾花6000购买过，不同的单细胞培训班，买的是一样的教程。现在免费送给你，别说话，开始学吧！
首先设置Seurat对象
我们将从 分析10X 外周血单核细胞 （PBMC） 数据集。原始数据可以在这里[1]找到。
读取数据。该函数从 10X 读取，返回独特的分子识别 （UMI） 计数矩阵。此矩阵中的值表示每个功能（即基因;行）在每个细胞（列）中检测到的分子数量。
我们接下来使用计数矩阵创建一个对象。
library(dplyr)
library(Seurat)
library(patchwork)
# Load the PBMC dataset
pbmc.data <- Read10X(data.dir = "../data/pbmc3k/filtered_gene_bc_matrices/hg19/")
# Initialize the Seurat object with the raw (non-normalized data).
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
pbmc
## An object of class Seurat
## 13714 features across 2700 samples within 1 assay
## Active assay: RNA (13714 features, 0 variable features)
标准预处理工作流程
以下步骤包括 Seurat 中 scRNA-seq 数据的标准预处理工作流程。包括基于 QC 指标的过滤、数据标准化和归一化，以及检测高变异基因的功能。
QC 和选择细胞以供进一步分析
Seurat 允许您轻松地探索 QC 指标，并根据任何用户定义的标准过滤细胞。常用[2]的一些 QC 指标包括
- 在每个细胞中检测到的基因数量。
- 低质量细胞或空液滴通常只有很少的基因
- 细胞双倍或多胞可能表现出异常高的基因计数
- 同样，细胞内检测到的分子总数（与独特的基因密切相关）
- 读取该细胞中的线粒体基因组的百分比
- 低质量/死细胞经常表现出广泛的线粒体污染
- 我们使用PercentageFeatureSet()[3]函数计算线粒体 QC 指标，该函数计算来自一组功能的计数百分比
# The [[ operator can add columns to object metadata. This is a great place to stash QC stats
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
在下面的示例中，我们可视化 QC 指标，并使用这些指标来过滤细胞。
- 过滤具有UMI计数超过 2，500 或小于 200的细胞
- 过滤具有>5%线粒体的细胞
# Visualize QC metrics as a violin plot
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
# FeatureScatter is typically used to visualize feature-feature relationships, but can be used
# for anything calculated by the object, i.e. columns in object metadata, PC scores etc.
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
标准化数据
从数据集中删除不需要的细胞后，下一步是使数据标准化。默认情况下，我们采用全局标准化。
pbmc <- NormalizeData(pbmc)
高变异基因的选择
接下来，我们计算数据集中显示高变异的特征子集（即，它们在某些细胞中表达强烈，在另一些单元格中表达得很低）。在下游分析中关注这些基因有助于突出单细胞数据集中的生物信号。
默认情况下，我们使用每个数据集的 2，000 个基因。这些将用于下游分析，如 PCA。
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
# Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(pbmc), 10)
# plot variable features with and without labels
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2
归一化数据
接下来，我们应用线性转换（"归一化"）ScaleData()[4]继续处理数据集。目的是
- 改变每个基因的表达，使细胞的平均表达为0
- 缩放每个基因的表达，使细胞之间的方差为 1
- 这一步骤在下游分析中具有同等的权重，因此高度表达的基因不会占主导地位
- 结果存储在
pbmc[["RNA"]]@scale.data
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
PCA分析
接下来，我们将在缩放数据上执行PCA。默认情况下，只有先前确定的可变功能用作输入，但如果您希望选择不同的子集，则可以使用参数进行定义。
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
Seurat 提供了几个有用的方法来可视化定义 PCA 的单元格和功能，包括 ，[DimPlot()]，[DimHeatmap()]等
# Examine and visualize PCA results a few different ways
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)
## PC_ 1
## Positive: CST3, TYROBP, LST1, AIF1, FTL
## Negative: MALAT1, LTB, IL32, IL7R, CD2
## PC_ 2
## Positive: CD79A, MS4A1, TCL1A, HLA-DQA1, HLA-DQB1
## Negative: NKG7, PRF1, CST7, GZMB, GZMA
## PC_ 3
## Positive: HLA-DQA1, CD79A, CD79B, HLA-DQB1, HLA-DPB1
## Negative: PPBP, PF4, SDPR, SPARC, GNG11
## PC_ 4
## Positive: HLA-DQA1, CD79B, CD79A, MS4A1, HLA-DQB1
## Negative: VIM, IL7R, S100A6, IL32, S100A8
## PC_ 5
## Positive: GZMB, NKG7, S100A8, FGFBP2, GNLY
## Negative: LTB, IL7R, CKB, VIM, MS4A7
VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
DimPlot(pbmc, reduction = "pca")
特别是，它允许在数据集中轻松探索异质性的主要来源，并且在尝试决定将哪些 PC 用于进一步下游分析
DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)
DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)
确定数据集的"主成分个数"
Seurat 根据 PCA 分数对单元单元进行聚类，每台 PC 基本上代表一个"元结构"，该"元结构"将信息组合在相关功能集中。我们随机排列数据的子集（默认情况下为 1%）并重新运行 PCA，构建功能分数的"空分布"，并重复此过程。我们确定"重要"PC。
# NOTE: This process can take a long time for big datasets, comment out for expediency. More
# approximate techniques such as those implemented in ElbowPlot() can be used to reduce
# computation time
pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
JackStrawPlot()功能提供了一个可视化工具，用于将每个 PC 的 p 值分布与统一分布（虚线）进行比较。"重要"PC 将显示在虚线上方。
JackStrawPlot(pbmc, dims = 1:15)
另一种启发式方法生成"肘部图"：[ElbowPlot()]根据每个（函数）解释的方差百分比对原则组件进行排名。在此示例中，我们可以观察到 PC9-10 周围的"肘部"，表明大多数真实信号在前 10 个 PC 中被捕获。
ElbowPlot(pbmc)
我们在这里选择了 10 个，但这不是固定的。
细胞聚类
Seurat 采用基于图形的聚类方法，简言之，这些方法将细胞嵌入到图形结构中，例如 K 最近的邻邻 （KNN） 图，在具有类似特征表达模式的单元之间绘制边缘，然后尝试将此图划分为高度互连的"集团"或"社区"。
与表象一样，我们首先根据 PCA 空间中的欧几里德距离构建 KNN 图，并根据当地社区的共享重叠（Jaccard 相似性）优化任意两个细胞之间的边缘权重。接下来将 Louvain 算法（默认值）或 SLM 等模块化优化技术应用于迭代组细胞，以优化标准模块化功能。该函数实现此过程，并包含一个分辨率参数，该参数设置下游聚类的"数量"，增加值导致更多群集。我们发现，将此参数设置在 0.4-1.2 之间通常会为大约 3K 细胞的单细胞数据集提供良好的结果。对于较大的数据集，最佳分辨率通常会增加。
pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 2638
## Number of edges: 95965
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8723
## Number of communities: 9
## Elapsed time: 0 seconds
# Look at cluster IDs of the first 5 cells
head(Idents(pbmc), 5)
## AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1 AAACCGTGCTTCCG-1
## 2 3 2 1
## AAACCGTGTATGCG-1
## 6
## Levels: 0 1 2 3 4 5 6 7 8
运行非线性降维（UMAP/tSNE）
Seurat 提供了多种非线性降维技术，如 tSNE 和 UMAP，以可视化和探索这些数据集。。
# If you haven't installed UMAP, you can do so via reticulate::py_install(packages =
# 'umap-learn')
pbmc <- RunUMAP(pbmc, dims = 1:10)
# note that you can set `label = TRUE` or use the LabelClusters function to help label
# individual clusters
DimPlot(pbmc, reduction = "umap")
此时，您可以保存对象，以便轻松加载，而无需重新运行上面执行的计算密集级步骤，或轻松与协作者共享。
saveRDS(pbmc, file = "../output/pbmc_tutorial.rds")
查找不同表达的marker
默认情况下，FindAllMarkers()与所有其他细胞相比，可识别单个群集的marker。但您也可以测试组群相互对立，或针对所有亚群。
# find all markers of cluster 2
cluster2.markers <- FindMarkers(pbmc, ident.1 = 2, min.pct = 0.25)
head(cluster2.markers, n = 5)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## IL32 2.593535e-91 1.2154360 0.949 0.466 3.556774e-87
## LTB 7.994465e-87 1.2828597 0.981 0.644 1.096361e-82
## CD3D 3.922451e-70 0.9359210 0.922 0.433 5.379250e-66
## IL7R 1.130870e-66 1.1776027 0.748 0.327 1.550876e-62
## LDHB 4.082189e-65 0.8837324 0.953 0.614 5.598314e-61
# find all markers distinguishing cluster 5 from clusters 0 and 3
cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)
head(cluster5.markers, n = 5)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## FCGR3A 2.150929e-209 4.267579 0.975 0.039 2.949784e-205
## IFITM3 6.103366e-199 3.877105 0.975 0.048 8.370156e-195
## CFD 8.891428e-198 3.411039 0.938 0.037 1.219370e-193
## CD68 2.374425e-194 3.014535 0.926 0.035 3.256286e-190
## RP11-290F20.3 9.308287e-191 2.722684 0.840 0.016 1.276538e-186
# find markers for every cluster compared to all remaining cells, report only the positive ones
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)
## # A tibble: 18 x 7
## # Groups: cluster [9]
## p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene
## <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <chr>
## 1 1.74e-109 1.07 0.897 0.593 2.39e-105 0 LDHB
## 2 1.17e- 83 1.33 0.435 0.108 1.60e- 79 0 CCR7
## 3 0\. 5.57 0.996 0.215 0\. 1 S100A9
## 4 0\. 5.48 0.975 0.121 0\. 1 S100A8
## 5 7.99e- 87 1.28 0.981 0.644 1.10e- 82 2 LTB
## 6 2.61e- 59 1.24 0.424 0.111 3.58e- 55 2 AQP3
## 7 0\. 4.31 0.936 0.041 0\. 3 CD79A
## 8 9.48e-271 3.59 0.622 0.022 1.30e-266 3 TCL1A
## 9 1.17e-178 2.97 0.957 0.241 1.60e-174 4 CCL5
## 10 4.93e-169 3.01 0.595 0.056 6.76e-165 4 GZMK
## 11 3.51e-184 3.31 0.975 0.134 4.82e-180 5 FCGR3A
## 12 2.03e-125 3.09 1 0.315 2.78e-121 5 LST1
## 13 1.05e-265 4.89 0.986 0.071 1.44e-261 6 GZMB
## 14 6.82e-175 4.92 0.958 0.135 9.36e-171 6 GNLY
## 15 1.48e-220 3.87 0.812 0.011 2.03e-216 7 FCER1A
## 16 1.67e- 21 2.87 1 0.513 2.28e- 17 7 HLA-DPB1
## 17 7.73e-200 7.24 1 0.01 1.06e-195 8 PF4
## 18 3.68e-110 8.58 1 0.024 5.05e-106 8 PPBP
几个可视化marker表达的工具。
VlnPlot(pbmc, features = c("MS4A1", "CD79A"))
# you can plot raw counts as well
VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = TRUE)
FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP",
"CD8A"))
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)
DoHeatmap(pbmc, features = top10$gene) + NoLegend()
亚群重命名
我们可以使用marker来轻松地将聚类与已知的单元类型进行匹配：
Cluster ID Markers Cell Type
- 0 IL7R, CCR7 Naive CD4+ T
- 1 CD14, LYZ CD14+ Mono
- 2 IL7R, S100A4 Memory CD4+
- 3 MS4A1 B
- 4 CD8A CD8+ T
- 5 FCGR3A, MS4A7 FCGR3A+ Mono
- 6 GNLY, NKG7 NK
- 7 FCER1A, CST3 DC
- 8 PPBP Platelet
new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "B", "CD8 T", "FCGR3A+ Mono",
"NK", "DC", "Platelet")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
最后保存数据。
saveRDS(pbmc, file = "../output/pbmc3k_final.rds")
参考资料
[1]
在这里: https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz
[2]
常用: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4758103/
[3]
PercentageFeatureSet(): https://satijalab.org/seurat/reference/PercentageFeatureSet.html
[4]
ScaleData(): https://satijalab.org/seurat/reference/ScaleData.html
往期回顾
OSCA单细胞数据分析笔记10—Marker gene detection
NC单细胞文章复现(七)：Gene expression signatures(2)
如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程
