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

引言

本系列开启 R 中单细胞RNA-seq数据分析教程，持续更新，欢迎关注，转发！

• 利用注释好的参考数据集辅助新数据分析

随着全球范围内单细胞RNA测序（scRNA-seq）数据的不断增多，特别是在人类细胞图谱（HCA）项目的推动下，大量注释详尽的图谱级scRNA-seq数据集已公开可用。因此，在分析新的相关数据集时，若不利用这些资源来辅助分析，尤其是帮助注释，将是一种资源浪费。这与之前的情况不同，之前处理的是多个地位平等的数据集整合。

在这里，有一个已经注释好的参考数据集和一个待分析的查询数据集。当然，仍可采用前文提到的相同方法，将参考数据集与查询数据集整合，随后进行参考-查询联合分析或以查询为中心的分析。但更直接的方法是进行数据转移或投影分析，即固定参考数据集，并尝试将查询数据集中的细胞或细胞群体与参考数据集中的对应部分进行匹配。

在本教程的这一部分，将介绍两到三种策略，说明如何利用参考数据来辅助注释查询数据。将使用前文提到的DS1作为查询数据，并从同一篇论文中选取一个注释好的数据子集作为参考数据。参考数据可通过此链接获取。在参考数据的元数据框架中，有一列名为“celltype”，记录了每个细胞的注释类型。

加载数据

首先导入Seurat工具及所需数据，包括已保存的DS1的Seurat对象和新下载的参考数据集的Seurat对象。

library(Seurat)
seurat_DS1 <- readRDS("DS1/seurat_obj_all.rds")
seurat_ref <- readRDS("ref_seurat_obj.rds")

来看看参考数据集是什么样子的

library(patchwork)
plot1 <- UMAPPlot(seurat_ref, group.by="branch")
plot2 <- UMAPPlot(seurat_ref, group.by="celltype")
plot3 <- FeaturePlot(seurat_ref, c("SOX2","DCX","FOXG1","EMX1","DLX2","LHX9"),
                     ncol=3, pt.size = 0.1)
((plot1 / plot2) | plot3) + plot_layout(width = c(1,3))

可以看到，参考数据集已经经过了完善的注释，其中包含了代表不同大脑区域和神经元亚型的多种细胞类型。

基于细胞簇的转录组相似性分析

可以将查询数据集中每个细胞群体的转录组特征与参考数据集中不同细胞类型的转录组特征进行比较。通过这种比较，查询数据中的细胞簇可以被归类为参考数据集中转录组特征最相似的细胞类型。为了实现这一目标，需要首先明确以下两点：

1. 选择用于计算转录组相似性的基因：
   一个直接的方法是使用参考数据集中的高变基因。还可以将这些高变基因与查询数据集中的可变基因取交集。此外，也可以先识别参考数据集中每种细胞类型的标记基因，然后将这些标记基因合并为一个代表转录组特征的基因列表。
2. 定义转录组特征之间的相似性：
   常用的方法是计算基因间的相关系数。常见的相关系数包括 Pearson 相关系数和 Spearman 相关系数。Pearson 相关系数主要用于衡量两个向量之间的线性关系，而 Spearman 相关系数则是通过将向量排序后计算 Pearson 相关系数，从而能够评估两个向量之间的非线性单调关系。Spearman 相关系数通常对少量异常值更具鲁棒性，因此能提供更稳健的估计结果。然而，排序操作通常耗时且计算资源消耗较大，尤其是在处理高通量数据时，这使得 Spearman 相关系数的计算速度远慢于 Pearson 相关系数。

在接下来的示例中，将使用参考数据集和查询数据集的高变基因的交集，通过计算 Spearman 相关系数来表示转录组相似性。

首先，需要分别计算参考数据集中每种注释细胞类型和查询数据集中每个细胞簇的平均转录组特征。

avg_expr_ref <- sapply(sort(unique(seurat_ref$celltype)), function(ct) rowMeans(seurat_ref@assays$RNA@data[,which(seurat_ref$celltype == ct)] ))
avg_expr_ds1 <- sapply(levels([e-mail]), function(ct) rowMeans(seurat_ds1@assays$RNA@data[,which([e-mail] == ct)]))

然后，找出能代表转录组的基因，再算出参考数据中各类细胞和查询数据中各细胞簇，这些基因平均表达水平之间的斯皮尔曼相关性，也就是两两之间的相关程度。

genes2cor <- intersect(VariableFeatures(seurat_ref), rownames(seurat_ds1))
corr2ref_cl <- cor(avg_expr_ds1[genes2cor,], avg_expr_ref[genes2cor,], method="spearman")

现在可以使用热图来可视化相关矩阵。

library(gplots)
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))

通过查看热图，就能依据转录组和参考数据里已标注细胞类型的相似度，来评估之前给这个数据集做的标注是否合理。比如查询数据里的不少细胞簇，像被标注为“Dorsal telen. IP”的那个细胞簇，就和参考数据中的“Cortical IP”细胞类型，有着最高的相似度。

本文由mdnice多平台发布