R中单细胞RNA-seq分析教程 (12)

引言

本系列开启 R 中scRNA-seq数据分析教程，持续更新，欢迎关注，转发！想要获取更多教程内容或者生信分析服务可以添加文末的学习交流群或客服QQ:941844452。

细胞层面的转录组相似性

前面文章提到的第一种方法试图把两个数据集中的细胞群或细胞类型关联起来。这种方法虽然简单，但有个明显的缺点：两个数据集里的细胞群或类型可能分辨率不同，没法直接比较。这在动态系统里尤其关键，比如研究发育或再生时，细胞状态是连续变化的，而聚类分析可能会在不同数据集里以不同方式切分这种连续性。碰到这种情况，一个替代办法是计算细胞与不同参考细胞类型的转录组相似性，但不是针对每个查询细胞群，而是对每个单独的查询细胞进行分析。

同样地，用两个数据集中高度可变基因的交集作为转录组的特征。至于用哪种相关性分析，因为单细胞转录组数据非常稀疏， Spearman 相关性通常效果更好。不过，像前面提到的，算 Spearman 相关性需要先对数据排序。给几千个基因在几千个细胞里的表达排序，不仅耗时间，还会生成一个巨大的密集矩阵，占内存很严重，细胞数量太多时甚至可能超过 R 环境的承受能力。所以，不能直接用 R 的基础函数 cor 来算。得找个更聪明的方法，在排序时保持矩阵稀疏，然后用专门处理稀疏矩阵的包，基于排序后的稀疏矩阵算 Pearson 相关性。

首先，对于参考数据集里细胞类型的平均表达profile，因为数据不稀疏，而且条目数（也就是细胞类型数量）不会太多，可以直接用基础的 rank 函数来排序。

ranked_expr_ref <- apply(avg_expr_ref[genes2cor,],2,rank)

接下来，会讲两种快速给稀疏矩阵排序的方法。第一种是用 presto 包里的 rank_matrix 函数。这个函数用 C 语言写了个快速排序算法，所以速度特别快。

library(presto)
ranked_expr_ds1 <- rank_matrix(seurat_DS1@assays$RNA@data[genes2cor,])$X_ranked

另外，也可以自己动手写一个稀疏排序算法。

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,])

最后，要快速算出两个稀疏矩阵之间，或者一个稀疏矩阵和一个密集矩阵之间的 Pearson 相关性，建议用 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

来对比一下之前的注释和这次的投影注释。

plot1 <- UMAPPlot(seurat_DS1, label=T)
plot2 <- UMAPPlot(seurat_DS1, group.by="celltype_maxcor", label=T)
plot1 | plot2

动动您发财的小手点个赞吧！欢迎转发！

本文由mdnice多平台发布