单细胞与空间转录组的跨模态标签转移

本文代码来自Seurat官方网站。介绍单样本空转数据的标签转移。

rm(list = ls())
library(tidyverse)
library(Seurat)
library(SeuratData)
cortex = readRDS("cortex.rds")

这是小鼠皮层细胞的高质量参考数据，来自 Allen 研究所。

下载自 https://www.dropbox.com/s/cuowvm4vrf65pvq/allen_cortex.rds?dl=1

allen_reference <- readRDS("allen_cortex.rds")
dim(allen_reference)

## [1] 34617 14249

设置 ncells=3000 会将整个数据集归一化，但仅在 3000 个细胞上学习噪声模型。 这能显著加快 SCTransform 的速度，且性能无损失。是官网推荐的用法。

allen_reference <- SCTransform(allen_reference, ncells = 3000, verbose = FALSE) %>%
                      RunPCA(verbose = FALSE) %>%
                      RunUMAP(dims = 1:30)

参考数据是已经注释好的Seurat对象，注释的细胞类型存储在meta.data的subclass列，画个umap图展示它。

DimPlot(allen_reference, group.by = "subclass", label = TRUE)+NoLegend()+coord_fixed()

简单罗列一下这些细胞类型咯，不研究这个领域的话，只要分清楚这些单词是细胞的名字，不是基因名字就可以啦！

兴奋性神经元	L2/3 IT	第2/3层端脑内投射神经元
兴奋性神经元	L4	第4层颗粒神经元
兴奋性神经元	L5 IT	第5层端脑内投射神经元
兴奋性神经元	L5 PT	第5层锥体束投射神经元
兴奋性神经元	L6 CT	第6层皮层丘脑投射神经元
兴奋性神经元	L6 IT	第6层端脑内投射神经元
兴奋性神经元	L6b	第6b层神经元
兴奋性神经元	NP	近距离投射神经元
抑制性神经元	Lamp5	Lamp5 中间神经元
抑制性神经元	Meis2	Meis2 中间神经元
抑制性神经元	Pvalb	小清蛋白中间神经元
抑制性神经元	Serpinf1	Serpinf1 中间神经元
抑制性神经元	Sncg	Sncg 中间神经元
抑制性神经元	Sst	生长抑素中间神经元
抑制性神经元	Vip	VIP 中间神经元
发育相关神经元	CR	Cajal-Retzius 细胞
非神经元细胞	Astro	星形胶质细胞
非神经元细胞	Endo	血管内皮细胞
非神经元细胞	Macrophage	巨噬细胞/小胶质细胞
非神经元细胞	Oligo	少突胶质细胞
非神经元细胞	Peri	周细胞
非神经元细胞	SMC	平滑肌细胞
非神经元细胞	VLMC	血管及软脑膜细胞

寻找锚点，进行标签转移。这里的 cortex 已经完成了 SCTransform 和 PCA，因此可以直接用于寻找锚点和标签转移。

anchors <- FindTransferAnchors(reference = allen_reference,
                               query = cortex,
                               normalization.method = "SCT",
                               reference.reduction = "pca",
                               dims = 1:30)

predictions.assay <- TransferData(anchorset = anchors,
                                  refdata = allen_reference$subclass, 
                                  prediction.assay = TRUE,
                                  weight.reduction = cortex[["pca"]], 
                                  dims = 1:30)

此时我们得到的predictions.assay是每个spot是每一种细胞的预测分数。把他插入空转Seurat对象作为一个组成部分。

GetAssayData(predictions.assay, layer = "data")[1:4, 1:4]

##       AAACAGAGCGACTCCT-1 AAACCGGGTAGGTACC-1 AAACCGTTCGTCCAGG-1
## Vip                    0                  0                  0
## Lamp5                  0                  0                  0
## Sst                    0                  0                  0
## Sncg                   0                  0                  0
##       AAACTCGTGATATAAG-1
## Vip                    0
## Lamp5                  0
## Sst                    0
## Sncg                   0

cortex[["predictions"]] <- predictions.assay
DefaultAssay(cortex) <- "predictions"

画图查看其中两种细胞的预测分数。

SpatialFeaturePlot(cortex, features = c("L2/3 IT", "L4"), ncol = 2, crop = TRUE) &
  theme(plot.margin = unit(c(1, 1, 1, 1), "mm"),
        legend.text = element_text(size = 8))

寻找(预测分数)有明显空间分布模式的细胞类型。

cortex <- FindSpatiallyVariableFeatures(cortex, 
                                        assay = "predictions", 
                                        selection.method = "moransi",
                                        features = rownames(cortex), 
                                        layer = "data")

展示 Moran’s I 最高的 4 种细胞类型，也就是预测分数空间自相关性最强、空间分布模式最明显的细胞类型。（图上画的是细胞类型预测分数，不是基因表达量，也不是真实细胞数量或者丰度。）

top.clusters <- head(SpatiallyVariableFeatures(cortex, method = "moransi"), 4)

SpatialFeaturePlot(object = cortex, features = top.clusters, ncol = 2) & 
  theme(plot.margin = unit(c(1, 1, 1, 1), "mm"),
        legend.text = element_text(size = 8))

最后，展示各种细胞类型预测分数的空间分布，可以与已知的小鼠皮层分层结构进行比较，检查是否相符。

features_of_interest <- c("Astro", "L2/3 IT", "L4", "L5 PT", "L5 IT", "L6 CT", "L6 IT", "L6b", "Oligo")

SpatialFeaturePlot(cortex, 
                   features = features_of_interest,
                   pt.size.factor = 1.6, ncol = 3, 
                   crop = TRUE,alpha = c(0.1, 1)) & 
  theme(plot.margin = unit(c(1, 1, 1, 1), "mm"))