查看: 727|回复: 7

[流程分析] NBIS系列单细胞转录组数据分析实战(四):细胞聚类分析

[复制链接]

版主

Rank: 10Rank: 10Rank: 10

主题
29
注册时间
2016.5.11
在线时间
558 小时

活跃会员突出贡献


发表于 2021.4.6 11:29:46 | 显示全部楼层 |阅读模式
本文首发于微信公众号“bioinfomics”:https://mp.weixin.qq.com/s/WW0_80OM6K_iTI2totH-HQ

第四节:细胞聚类分析
在本节教程中,我们将基于批次矫正后的整合数据集进行细胞聚类分析,我们使用PCA线性降维的结果分别执行k-最近邻图聚类,层次聚类和k-均值聚类。

加载所需的R包和数据集
[AppleScript] 纯文本查看 复制代码
if (!require(clustree)) {
    install.packages("clustree", dependencies = FALSE)
}
## Loading required package: clustree
## Loading required package: ggraph

suppressPackageStartupMessages({
    library(Seurat)
    library(cowplot)
    library(ggplot2)
    library(pheatmap)
    library(rafalib)
    library(clustree)
})

alldata <- readRDS("data/results/covid_qc_dr_int.rds")


执行k-最近邻图聚类
在执行图聚类的过程中主要包括以下3个步骤:
  • Build a kNN graph from the data
  • Prune spurious connections from kNN graph (optional step). This is a SNN graph.
  • Find groups of cells that maximizes the connections within the group compared other groups.

构建kNN/SNN图
执行图聚类的第一步是构建一个kNN图,我们使用PCA降维的前N个PC用于计算。
我们可以使用Seurat包中的FindNeighbors函数计算构建KNN和SNN图。
[AppleScript] 纯文本查看 复制代码
# check that CCA is still the active assay
[url=mailto:alldata@active.assay]alldata@active.assay[/url]
## [1] "CCA"

# 使用FindNeighbors函数构建SNN图
alldata <- FindNeighbors(alldata, dims = 1:30, k.param = 60, prune.SNN = 1/15)
## Computing nearest neighbor graph
## Computing SNN

# check the names for graphs in the object.
names(alldata@graphs)
## [1] "CCA_nn"  "CCA_snn"

我们可以看一下kNN图,它是一个连接矩阵,其中不同细胞之间的每个连接都表示为1个s,这称之为未加权图(Seurat中的默认值)。但是,某些细胞之间的连接可能比其他细胞的更重要,在这种情况下,图的尺度会从0到最大距离。通常,距离越小,两点越接近,它们之间的连接也越牢固,这称之为加权图。加权图和未加权图均适用于图聚类,但是对于大型数据集(>100k细胞),使用非加权图在聚类上的速度会更快。
[AppleScript] 纯文本查看 复制代码
# check that CCA is still the active assay
[url=mailto:alldata@active.assay]alldata@active.assay[/url]
## [1] "CCA"

# 使用FindNeighbors函数构建SNN图
alldata <- FindNeighbors(alldata, dims = 1:30, k.param = 60, prune.SNN = 1/15)
## Computing nearest neighbor graph
## Computing SNN

# check the names for graphs in the object.
names(alldata@graphs)
## [1] "CCA_nn"  "CCA_snn"

基于SNN图进行细胞聚类
在构建好SNN图后,我们可以基于其执行图聚类。选用不同的分辨率(resolution)进行细胞聚类,分辨率越大,聚类出来的细胞簇数越多。
在Seurat中,我们使用FindClusters函数进行细胞聚类,默认情况下(algorithm = 1),该函数将使用“ Louvain”算法进行基于图的聚类。要使用leiden算法,我们需要将其设置为algorithm = 4。
[AppleScript] 纯文本查看 复制代码
# Clustering with louvain (algorithm 1)
for (res in c(0.1, 0.25, 0.5, 1, 1.5, 2)) {
    alldata <- FindClusters(alldata, graph.name = "CCA_snn", resolution = res, algorithm = 1)
}

# each time you run clustering, the data is stored in meta data columns:
# seurat_clusters - lastest results only CCA_snn_res.XX - for each different
# resolution you test.

plot_grid(ncol = 3, DimPlot(alldata, reduction = "umap", group.by = "CCA_snn_res.0.5") + ggtitle("louvain_0.5"), 
                    DimPlot(alldata, reduction = "umap", group.by = "CCA_snn_res.1") + ggtitle("louvain_1"), 
                    DimPlot(alldata, reduction = "umap", group.by = "CCA_snn_res.2") + ggtitle("louvain_2"))

现在,我们可以使用clustree包来可视化不同分辨率下细胞在聚类群之间的分配。
[AppleScript] 纯文本查看 复制代码
# install.packages('clustree')
suppressPackageStartupMessages(library(clustree))

clustree([url=mailto:alldata@meta.data]alldata@meta.data[/url], prefix = "CCA_snn_res.")

K均值聚类
K-means是一种常用的聚类算法,已在许多应用领域中使用。在R中,可以通过kmeans函数进行调用。通常,它应用于表达数据的降维表示(由于低维距离的可解释性,因此通常用于PCA)。

我们需要预先设定聚类群的数量。由于聚类的结果取决于群集中心的初始化,因此通常建议使用多个启动配置(通过nstart参数)运行K-means。
[AppleScript] 纯文本查看 复制代码
for (k in c(5, 7, 10, 12, 15, 17, 20)) {
    [url=mailto:alldata@meta.data]alldata@meta.data[/url][, paste0("kmeans_", k)] <- kmeans(x = alldata@reductions[["pca"]]@cell.embeddings,  centers = k, nstart = 100)$cluster
}

plot_grid(ncol = 3, DimPlot(alldata, reduction = "umap", group.by = "kmeans_5") +  ggtitle("kmeans_5"), 
                    DimPlot(alldata, reduction = "umap", group.by = "kmeans_10") +  ggtitle("kmeans_10"), 
                    DimPlot(alldata, reduction = "umap", group.by = "kmeans_15") +  ggtitle("kmeans_15"))

使用clustree函数查看不同聚类群的结果
[AppleScript] 纯文本查看 复制代码
clustree([url=mailto:alldata@meta.data]alldata@meta.data[/url], prefix = "kmeans_")

层次聚类定义细胞之间的距离
基本的R stats包中包含一个dist函数,可以用于计算所有成对样本之间的距离。由于我们要计算样本之间的距离,而不是基因之间的距离,因此我们需要先对表达数据进行转置,然后再将其应用于dist函数中。dist函数中可用的距离计算方法有:“euclidean”, “maximum”, “manhattan”, “canberra”, “binary” or “minkowski”.
[AppleScript] 纯文本查看 复制代码
d <- dist(alldata@reductions[["pca"]]@cell.embeddings, method = "euclidean")

可以看到,dist函数不能实现correlation的方法。但是,我们可以创建自己的距离并将其转换为距离对象。我们首先可以使用cor函数计算样本之间的相关性。如您所知,相关性的范围是从-1到1的,其中1表示两个样本最接近,-1表示两个样本最远,0介于两者之间。但是,这在定义距离时会产生问题,因为距离0表示两个样本最接近,距离1表示两个样本最远,而距离-1没有意义。因此,我们需要将相关性转换为正尺度(又称adjacency):

将相关性转换为0-1比例后,我们可以简单地使用as.dist函数将其转换为距离对象。
[AppleScript] 纯文本查看 复制代码
# Compute sample correlations
# 计算细胞之间的相关性
sample_cor <- cor(Matrix::t(alldata@reductions[["pca"]]@cell.embeddings))

# Transform the scale from correlations
sample_cor <- (1 - sample_cor)/2

# Convert it to a distance object
d2 <- as.dist(sample_cor)

基于细胞之间的距离进行层次聚类
在计算出所有样本之间的距离之后,我们可以对其进行层次聚类。我们将使用hclust函数实现该功能,在该函数中,我们可以简单地使用上面创建的距离对象来运行它。可用的方法有:“ward.D”, “ward.D2”, “single”, “complete”, “average”, “mcquitty”, “median” or “centroid”。
[AppleScript] 纯文本查看 复制代码
# euclidean
h_euclidean <- hclust(d, method = "ward.D2")

# correlation
h_correlation <- hclust(d2, method = "ward.D2")

创建好分层聚类树后,下一步就是定义哪些样本属于特定簇。我们可以使用cutree函数根据特定k值切割聚类树,以定义聚类群。我们还可以定义簇的数量或确定高度。
[AppleScript] 纯文本查看 复制代码
#euclidean distance
alldata$hc_euclidean_5 <- cutree(h_euclidean,k = 5)
alldata$hc_euclidean_10 <- cutree(h_euclidean,k = 10)
alldata$hc_euclidean_15 <- cutree(h_euclidean,k = 15)

#correlation distance
alldata$hc_corelation_5 <- cutree(h_correlation,k = 5)
alldata$hc_corelation_10 <- cutree(h_correlation,k = 10)
alldata$hc_corelation_15 <- cutree(h_correlation,k = 15)

plot_grid(ncol = 3,
  DimPlot(alldata, reduction = "umap", group.by = "hc_euclidean_5")+ggtitle("hc_euc_5"),
  DimPlot(alldata, reduction = "umap", group.by = "hc_euclidean_10")+ggtitle("hc_euc_10"),
  DimPlot(alldata, reduction = "umap", group.by = "hc_euclidean_15")+ggtitle("hc_euc_15"),

  DimPlot(alldata, reduction = "umap", group.by = "hc_corelation_5")+ggtitle("hc_cor_5"),
  DimPlot(alldata, reduction = "umap", group.by = "hc_corelation_10")+ggtitle("hc_cor_10"),
  DimPlot(alldata, reduction = "umap", group.by = "hc_corelation_15")+ggtitle("hc_cor_15")
)

保存细胞聚类的结果
[AppleScript] 纯文本查看 复制代码
saveRDS(alldata, "data/results/covid_qc_dr_int_cl.rds")


本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
加油啦。
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
3
注册时间
2017.9.8
在线时间
28 小时

发表于 2021.4.7 08:42:27 | 显示全部楼层
加油,加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2020.11.21
在线时间
32 小时

发表于 2021.4.7 10:09:06 | 显示全部楼层
坚持就是胜利!
回复

使用道具 举报

草履虫

Rank: 2

主题
0
注册时间
2020.11.11
在线时间
1 小时

发表于 2021.4.7 10:39:23 | 显示全部楼层
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
3
注册时间
2017.9.8
在线时间
28 小时

发表于 2021.4.8 08:24:57 | 显示全部楼层
加油,加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2020.11.21
在线时间
32 小时

发表于 2021.4.8 10:38:24 | 显示全部楼层
坚持就是胜利!
回复

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2021.3.29
在线时间
3 小时

发表于 2021.4.9 11:17:57 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2021.3.29
在线时间
3 小时

发表于 2021.4.9 11:22:19 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

快速回复 返回顶部 返回列表