引言
本系列开启 R 中单细胞RNA-seq数据分析教程，持续更新，欢迎关注，转发！
细胞水平上的转录组相似性分析
第一种方法的目标是将两个数据集中的细胞簇或细胞类型进行关联。这种方法虽然简单，但存在一个明显的缺陷：两个数据集中的细胞簇或细胞类型可能没有以相同的分辨率定义，因此难以直接比较。这一问题在动态系统中尤为突出，例如发育或再生过程，这些系统中细胞状态是连续变化的，而聚类分析可能会因数据集的不同而以不同的方式分割这些连续状态。在这种情况下，一种替代方案是计算查询细胞与不同参考细胞类型的转录组相似性，但不是针对每个查询细胞簇进行计算，而是针对每个查询细胞单独计算。
我们采用两个数据集中高度可变基因的交集作为转录组特征。在选择相关性分析方法时，由于单细胞转录组数据高度稀疏，Spearman 相关性分析通常表现更好。然而，如前所述，计算 Spearman 相关性需要对数据进行排序。对数千个基因在数千个细胞中的表达进行排序不仅耗时，还会生成一个巨大的密集矩阵，这需要大量内存，有时甚至会超出 R 环境的处理能力，尤其是在细胞数量极多的情况下。因此，我们不能依赖 R 的基本函数 cor 来完成这一计算。我们需要一种更高效的方法来实现排序，同时保持矩阵的稀疏性，然后利用专门用于稀疏矩阵的 Pearson 相关性计算包，针对排序后的稀疏矩阵进行相关性分析。
对于参考数据集的细胞类型平均表达谱，由于其本身并不稀疏，且涉及的细胞类型数量通常不多，我们可以直接使用基本的 rank 函数进行处理。
ranked_expr_ref <- apply(avg_expr_ref[genes2cor,],2,rank)
AI 代码解读
接下来，我们将介绍两种稀疏矩阵快速排序的方法。第一种方法是利用 presto 包中的 rank_matrix 函数。这个包之前曾被用于标记物识别。该函数通过 C 语言实现了一种快速排序算法，因此运行速度极快。
library(presto)
ranked_expr_ds1 <- rank_matrix(seurat_DS1@assays$RNA@data[genes2cor,])$X_ranked
AI 代码解读
或者，我们也可以自行实现稀疏排序算法。你无需深究该函数的所有细节，只需进行复制粘贴操作。
rank_matrix <- function (mat)
{
if (is.matrix(mat) | is.data.frame(mat)) {
ranked_mat <- apply(mat, 2, rank)
}
else {
df_mat <- Matrix::summary(mat)
dfs_mat <- split(df_mat, df_mat$j)
df_mat_ranked <- do.call(rbind, lapply(dfs_mat, function(df) {
num_zeros <- nrow(mat) - nrow(df)
ranks_nonzero <- rank(df$x)
df$x <- ranks_nonzero + num_zeros - (1 + num_zeros)/2
return(df)
}))
ranked_mat <- sparseMatrix(i = df_mat_ranked$i, j = df_mat_ranked$j,
x = df_mat_ranked$x, dims = dim(mat), dimnames = dimnames(mat))
}
return(ranked_mat)
}
ranked_expr_ds1 <- rank_matrix(seurat_DS1@assays$RNA@data[genes2cor,])
AI 代码解读
最后，若需快速计算两个稀疏矩阵之间或一个稀疏矩阵与一个密集矩阵之间的皮尔逊相关性，建议使用 qlcMatrix 包中的 corSparse 函数。随后，如果查询数据集中某个细胞的转录组与某种细胞类型的相似性最高，就可以将该细胞归类为该参考细胞类型。
library(qlcMatrix)
corr2ref_cell <- corSparse(ranked_expr_ds1, ranked_expr_ref)
ct_maxcor <- colnames(avg_expr_ref)[apply(corr2ref_cell, 1, which.max)]
seurat_DS1$celltype_maxcor <- ct_maxcor
AI 代码解读
我们来对比一下之前的注释结果和这次的预测注释结果。
plot1 <- UMAPPlot(seurat_DS1, label=T)
plot2 <- UMAPPlot(seurat_DS1, group.by="celltype_maxcor", label=T)
plot1 | plot2
AI 代码解读
虽然结果并不完美，但看起来也还算不错。
我们可以通过计算查询细胞簇中细胞与不同参考细胞类型的缩放相似性的平均值，来总结细胞水平的相似性，并将其以热图的形式可视化，类似于前面展示的那样。
corr2ref_scaled <- scale(t(corr2ref_cell))
corr2ref_sum2cl <- t(sapply(levels(seurat_DS1@active.ident), function(cl)
rowMeans(corr2ref_scaled[,which(seurat_DS1@active.ident == cl)]) ))
heatmap.2(corr2ref_cl, scale="none", trace="none", key=F, keysize=0.5, margins=c(15,17),
labRow = colnames(avg_expr_ds1), labCol = colnames(avg_expr_ref), cexRow=0.8, cexCol=0.8,
col=colorRampPalette(rev(c("#b2182b","#d6604d","#f4a582","#fddbc7","#f7f7f7","#d1e5f0","#92c5de","#4393c3","#2166ac")))(30))
AI 代码解读
