1. 加载包并读入数据
library(dplyr)
library(Seurat)
library(patchwork)
#读入数据
pbmc.data <- Read10X(data.dir = "D:R program /filtered_gene_bc_matrices/hg19/")
#创建seurat对象-未经标准化的
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
2. 查看数据信息
pbmc
# An object of class Seurat
# 13714 features across 2700 samples within 1 assay
# Active assay: RNA (13714 features, 0 variable features)
pbmc.data[c("CD3D", "TCL1A", "MS4A1"), 1:30]#查看"CD3D", "TCL1A", "MS4A1"这三行的1到30列
# 3 x 30 sparse Matrix（稀疏矩阵） of class "dgCMatrix"
#
# CD3D 4 . 10 . . 1 2 3 1 . . 2 7 1 . . 1 3 . 2 3 . . . . . 3 4 1 5
# TCL1A . . . . . . . . 1 . . . . . . . . . . . . 1 . . . . . . . .
# MS4A1 . 6 . . . . . . 1 1 1 . . . . . . . . . 36 1 2 . . 2 . . . .
矩阵中的.值表示 0（未检测到分子）。由于 scRNA-seq 矩阵中的大多数值为 0，Seurat 尽可能使用稀疏矩阵表示。这会显着节省 Drop-seq/inDrop/10x 数据的内存和速度。
3. 标准流程
- 读入数据
- 构建Seurat对象
- 预处理，QC
- 归一化
- 寻找特征基因
- 标准化
- PCA降维
- 分群
- 非线性降维UMAP和TSNE
- 可视化降维结果
pbmc.counts <- Read10X(data.dir = "D:R program /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")
3.3 质控
# 这个命令会在metadata中加入percent.mt这一列数据
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
-
在每个细胞中检测到的独特基因的数量。
- 低质量细胞或空液滴通常含有很少的基因
- 细胞双联体或多联体可能表现出异常高的基因计数
-
同样，细胞内检测到的分子总数（与独特基因密切相关）
-
映射到线粒体基因组的读数百分比
- 低质量/垂死的细胞通常表现出广泛的线粒体污染
- 我们使用函数计算线粒体 QC 指标，该
PercentageFeatureSet()函数计算源自一组特征的计数百分比 - 我们使用以 开头的所有基因
MT-的集合 作为一组线粒体基因
# Visualize QC metrics as a violin plot
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
# FeatureScatter可视化细胞与细胞之间的关系
# 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)
#这些内容是自定义的
3.4 归一化
从数据集中删除不需要的单元格后，下一步是规范化数据。默认情况下，我们采用全局尺度归一化方法“LogNormalize”，通过总表达式对每个单元格的特征表达式测量值进行归一化，将其乘以比例因子（默认为 10,000），并对结果进行对数转换。标准化值存储在pbmc[["RNA"]]@data.
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- NormalizeData(pbmc)#运行这个命令用的默认值，就是上面的LogNormalize和scale.factor = 10000
#这俩是一样的命令
3.5 寻找特征基因--高变基因
高细胞间变异的特征子集，它们在某些细胞中高表达，而在其他细胞中低表达，在下游分析中关注这些基因有助于突出单细胞数据集中的生物信号。
在这里详细描述了我们在 Seurat 中的过程，并通过直接建模单细胞数据中固有的均值-方差关系来改进以前的版本，并在FindVariableFeatures()函数中实现。默认情况下，我们为每个数据集返回 2,000 个特征。这些将用于下游分析，如 PCA。
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
# 10个高变基因
top10 <- head(VariableFeatures(pbmc), 10)
# 高变基因加标签-前10有标签
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2
3.6 标准化--缩放数据---主要用于热图展示
接下来，我们应用线性变换（“缩放”），这是 PCA 等降维技术之前的标准预处理步骤。缩放是 Seurat 工作流程中必不可少的一步，但仅限于将用作 PCA 输入的基因。ScaleData()功能：
-
移动每个基因的表达，使跨细胞的平均表达为 0
-
缩放每个基因的表达，使得跨细胞的方差为 1
- 此步骤在下游分析中给予同等权重，因此高表达基因不会占主导地位
-
其结果存储在
pbmc[["RNA"]]@scale.data
all.genes <- rownames(pbmc) #提取列作为所有的基因
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- ScaleData(pbmc)#这个跟上面的依据一个意思，都是默认的值2000
PCA 和聚类结果将不受影响。然而，Seurat 热图（如下图所示生成DoHeatmap()）需要热图中的基因进行缩放，以确保高表达基因不会主导热图。
3.7 线性降维
接下来我们对缩放后的数据执行 PCA。默认情况下，只有先前确定的变量特征用作输入，但如果选择不同的子集，可以使用features定义
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
Seurat 提供了几种有用的方法来可视化定义 PCA 的单元格和特征，包括VizDimReduction()、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")
可以看下批次效应
特别是DimHeatmap()允许轻松探索数据集中异质性的主要来源，并且决定要包括哪些 PC 以进行进一步的下游分析。单元格和特征都根据其 PCA 分数排序。设置cells为一个数字会在光谱两端绘制“极端”单元格，这会显着加快大型数据集的绘制速度。虽然显然是监督分析，但我们发现这是探索相关特征集的宝贵工具。
DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)
DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)
确定数据集的“维度”--选多不选少
顶级主成分代表了数据集的强大压缩
- 上面的DimHeatmap算一个
- 肘图算一个
- JackStrawPlot不太好看
在Macosko等人中，我们实施了受 JackStraw 程序启发的重采样测试。我们随机排列数据的一个子集（默认为 1%）并重新运行 PCA，构建特征分数的“零分布”，然后重复此过程。我们将“重要”PC 识别为具有低 p 值特征的大量丰富的 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(pbmc, dims = 1:15)
该JackStrawPlot()函数提供了一个可视化工具，用于将每个 PC 的 p 值分布与均匀分布（虚线）进行比较。“重要”PC 将显示具有低 p 值（虚线上方的实线）。在这种情况下，似乎在前 10-12 个 PC 之后显着性急剧下降。
ElbowPlot()另一种启发式方法生成“肘部图”：根据每个主成分（函数）解释的方差百分比对主成分进行排名。在此示例中，我们可以观察到 PC9-10 周围的“弯头”，这表明大部分真实信号是在前 10 个 PC 中捕获的。
ElbowPlot(pbmc)
pbmc <- FindNeighbors(pbmc, dims = 1:10)#KNN
pbmc <- FindClusters(pbmc, resolution = 0.5)#补充，Louvain 算法
# 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
3.8 非线性降维
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")
3.9 寻找差异表达的特征（簇生物标志物）
ident.1与所有其他细胞相比，它识别单个簇（在中指定）的阳性和阴性标记。FindAllMarkers()为所有集群自动执行此过程，也可以测试集群组之间的对比
min.pct参数要求在两组单元格中以最小百分比检测到一个特征，并且 thresh.test 参数要求一个特征在两组之间以一定数量差异表达（平均）。您可以将这两个设置为 0，但时间会急剧增加。
max.cells.per.ident可以设置用来加快时间。虽然通常会出现一部分数据损失，但速度可能会显着提高，并且差异表达最高的特征可能仍会上升到顶部。
# find all markers of cluster 2
cluster2.markers <- FindMarkers(pbmc, ident.1 = 2, min.pct = 0.25)
head(cluster2.markers, n = 5)
#ident 0，ident1,ident2分别显示阳性，阴性以及阳性和阴性的标记
# p_val avg_log2FC pct.1 pct.2 p_val_adj
# IL32 2.892340e-90 1.2013522 0.947 0.465 3.966555e-86
# LTB 1.060121e-86 1.2695776 0.981 0.643 1.453850e-82
# CD3D 8.794641e-71 0.9389621 0.922 0.432 1.206097e-66
# IL7R 3.516098e-68 1.1873213 0.750 0.326 4.821977e-64
# LDHB 1.642480e-67 0.8969774 0.954 0.614 2.252497e-63
# 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 8.246578e-205 4.261495 0.975 0.040 1.130936e-200
# IFITM3 1.677613e-195 3.879339 0.975 0.049 2.300678e-191
# CFD 2.401156e-193 3.405492 0.938 0.038 3.292945e-189
# CD68 2.900384e-191 3.020484 0.926 0.035 3.977587e-187
# RP11-290F20.3 2.513244e-186 2.720057 0.840 0.017 3.446663e-182
# 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) %>%
slice_max(n = 2, order_by = avg_log2FC)
# A tibble: 18 × 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 9.57e- 88 1.36 0.447 0.108 1.31e- 83 0 CCR7
# 2 3.75e-112 1.09 0.912 0.592 5.14e-108 0 LDHB
# 3 0 5.57 0.996 0.215 0 1 S100A9
# 4 0 5.48 0.975 0.121 0 1 S100A8
# 5 1.06e- 86 1.27 0.981 0.643 1.45e- 82 2 LTB
# 6 2.97e- 58 1.23 0.42 0.111 4.07e- 54 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 5.61e-202 3.10 0.983 0.234 7.70e-198 4 CCL5
## 10 7.25e-165 3.00 0.577 0.055 9.95e-161 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 3.13e-191 5.32 0.961 0.131 4.30e-187 6 GNLY
## 14 7.95e-269 4.83 0.961 0.068 1.09e-264 6 GZMB
## 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 1.92e-102 8.59 1 0.024 2.63e- 98 8 PPBP
## 18 9.25e-186 7.29 1 0.011 1.27e-181 8 PF4
Seurat 有几个差异表达检测方式，可以使用 test.use 参数设置（有关详细信息，请参阅DE vignette）。例如，ROC 测试返回任何单个标记的“分类能力”（范围从 0 - 随机到 1 - 完美）。
cluster0.markers <- FindMarkers(pbmc, ident.1 = 0, logfc.threshold = 0.25, test.use = "roc", only.pos = TRUE)
可视化标记表达工具
VlnPlot()（显示跨集群的表达概率分布）
FeaturePlot()（在 tSNE 或 PCA 图上可视化特征表达）是我们最常用的可视化。
也可以试着将RidgePlot()
CellScatter()
DotPlot()作为查看数据集的其他方法。
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"))
DoHeatmap()为给定的单元格和特征生成表达式热图。在这种情况下，我们为每个集群绘制前 20 个标记（或所有标记，如果少于 20 个）。
pbmc.markers %>%
group_by(cluster) %>%
top_n(n = 10, wt = avg_log2FC) -> top10
DoHeatmap(pbmc, features = top10$gene) + NoLegend()
3.10 细胞类型注释
| 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")
