Dimensional Reduction and Clustering
clp
06 June, 2020
加载前一步获取的R环境变量,
#load R libraries necessary
library(dplyr)
library(Seurat)
library(ggplot2)
load('filtered_gene_bc_matrices/hg19/01_pbmc3k_sctf.rd') #pbmc
降维
我们对缩放后的数据执行PCA聚类分析。默认情况下，只有先前确定的特征变量用作输入，但如果您希望选择不同的子集，则可以使用features参数进行自定义。
pbmc <- RunPCA(pbmc, verbose = FALSE)
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: FTL, LYZ, FTH1, CST3, S100A9
#> Negative: MALAT1, CCL5, NKG7, RPS27A, LTB
#> PC_ 2
#> Positive: NKG7, CCL5, GZMB, GNLY, GZMA
#> Negative: CD74, HLA-DRA, CD79A, HLA-DPB1, HLA-DQA1
#> PC_ 3
#> Positive: CD74, HLA-DRA, CD79A, NKG7, HLA-DPB1
#> Negative: S100A8, S100A9, JUNB, RPS12, IL7R
#> PC_ 4
#> Positive: S100A8, S100A9, LYZ, LGALS2, CD14
#> Negative: FCGR3A, LST1, FCER1G, AIF1, MS4A7
#> PC_ 5
#> Positive: GPX1, CCL5, PPBP, PF4, SDPR
#> Negative: FCGR3A, LST1, GNLY, AIF1, IFITM3
# Load PC1 and PC2 dimension vectors
VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
绘制2D PCA。请注意，PCA是根据SCT分析计算的。
DimPlot(pbmc, reduction = "pca")
具体地说，DimHeatmap允许轻松探索数据集中的主要异质性来源，并且在试图决定包括哪些PC以进行进一步的下游分析时非常有用。细胞和特征都根据它们的PCA得分进行排序。将cells设置为一个数字将绘制光谱两端的“极端”细胞，这大大加快了大型数据集的绘制速度。通过有监督的分析，我们发现这是探索相关特征值的一个有价值的工具。
DimHeatmap(pbmc, dims = 1:9, cells = 500, balanced = TRUE)
确定数据集的“维数”
为了克服scRNA-seq数据的任何单一特征中广泛的技术噪音，seurat根据细胞的PCA得分对细胞进行聚类，每个主成分PC实质上代表一个将相关特征集上的信息组合在一起的“元特征”。因此，主成分表示数据集的稳健压缩。但是，我们应该选择包含多少个组分？10个？20个？100个？
我们使用“Elbow plot”（拐点图）来确定适当数量的主成分。它根据每个主成分解释的方差百分比绘制主成分排名(ElbowPlot函数)。在本例中，我们可以观察到PC9-14周围的“拐点”，这表明大部分真实信号是在前13个PC中捕获的。
ElbowPlot(pbmc)
对于用户来说，识别数据集的真实维度可能是具有挑战性的/不确定的。检验标准差是一种启发式方法，但它是常用的，并且可以立即计算出来。我们可能有理由选择PC9-14之间的任何一个作为分界点
对细胞进行聚类
Seurat v3应用在Macosko et. al.发表的基于图片的聚类方法。重要的是，驱动聚类分析的距离度量(基于先前识别的PC)保持不变。然而，我们将蜂窝距离矩阵划分为类的方法有了显著改进。我们的方法受到了最近的手稿的很大启发，这些手稿将基于图的聚类方法应用于scRNA-seq数据CyTOF data (PhenoGraph, Levine et. al., Cell, 2015)。简而言之，这些方法将细胞嵌入到图结构中-例如K-最近邻近(KNN)图，在具有相似特征表达模式的单元之间绘制边，然后尝试将该图划分为高度互连的“quasi-cliques”或“communities”。
与PhenoGraph一样，我们首先基于PCA空间中的欧几里德距离构造KNN图，并基于任意两个单元在其局部邻域中共享的重叠(Jaccard相似度)来细化它们之间的边权重。此步骤使用FindNeighbors函数执行，并将先前定义的数据集维度(前13个PC)作为输入。
为了对细胞进行群集，我们接下来应用模块化优化技术，如Louvain算法(默认)或SLM (SLM, Blondel et. al., Journal of Statistical Mechanics)，将细胞迭代地分组在一起，目的是优化标准模块化功能。FindClusters函数实现此过程，并包含一个分辨率参数，该参数设置下游集群的“粒度”，值越大，集群数量越多。我们发现，将此参数设置在0.4-1.2之间对于大约3K个细胞的单细胞数据集通常会返回良好的结果。对于较大的数据集，最佳分辨率通常会增加。集群（clusters）可以通过Idents函数找到。
# Find K-nearest neigbors defined by the euclidean distance in PCA space
pbmc <- FindNeighbors(pbmc, dims = 1:13)
# Louvain algorithm to find communities such that the modularity is optimized
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: 94814
#>
#> Running Louvain algorithm...
#> Maximum modularity in 10 random starts: 0.8852
#> 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
#> 3 4 0 1
#> AAACCGTGTATGCG-1
#> 5
#> Levels: 0 1 2 3 4 5 6 7 8
DimPlot(pbmc, reduction = "pca", label = TRUE) + NoLegend()
#> Warning: Using `as.character()` on a quosure is deprecated as of rlang 0.3.0.
#> Please use `as_label()` or `as_name()` instead.
#> This warning is displayed once per session.
运行非线性降维(UMAP/tSNE)
Seurat提供了几种非线性降维技术，如tSNE和UMAP，以可视化和探索这些数据集。这些算法的目标是学习数据的底层流形，以便在低维空间中将相似的细胞放在一起。上述确定的基于图形的群集内的细胞应该共同定位在这些降维曲线图上。作为UMAP和tSNE的输入，我们建议使用相同的PC作为聚类分析的输入。
# In Seurat v3.1.2, RunUMAP calls `uwot` (umap implementation in R and it should be already installed) as default. If RunUMAP does not find R umap, then you can try Python umap via reticulate::py_install(packages = 'umap-learn')
pbmc <- RunUMAP(pbmc, dims = 1:13)
#> Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
#> To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
#> This message will be shown once per session
# note that you can set `label = TRUE` or use the LabelClusters function to help label
# individual clusters
DimPlot(pbmc, reduction = "umap", label = TRUE) + NoLegend()
让我们保存R变量，这样我们就可以继续工作了。
save(pbmc, file = "filtered_gene_bc_matrices/hg19/02_pbmc3k_cluster.rd",compress = TRUE)
本节重点
- PCA
- Clustering
- UMAP
