本文的代码来自Seurat官方网站的空转教程。主要的两个部分是：
- Harmony 批次校正：校正多个空间切片之间的技术差异，用校正后的低维空间进行邻居构建、聚类和 UMAP 展示。
- 单细胞标签转移：以带细胞类型注释的 scRNA-seq 数据为参考，为每个空间 spot 计算各细胞类型的预测分数。
需要注意，Visium spot 通常包含多个细胞。Seurat 标签转移给出的是参考细胞类型与每个 spot 的相似性分数，可以将最高分对应的细胞类型作为该 spot 的预测标签，不等同于细胞比例反卷积。
1. 加载包
rm(list = ls())
options(future.globals.maxSize = 10 * 1024^3)
set.seed(1234)
library(Seurat)
library(SeuratData)
library(harmony)
library(ggplot2)
library(patchwork)
数据来自10X官方网站，也可以从SeuratData获取。InstallData会安装R包，如果已安装，运行这句代码不会导致重复安装。
InstallData("stxBrain")
#install.packages("harmony")
2. 数据预处理
先对每个切片单独运行 SCTransform()，以便分别学习各切片的测序深度和技术噪声模型。samples 可以扩展为更多样本。
用merge把多个切片合并到一个Seurat对象中，才能放在一起做降维聚类分群。
#样本文件夹名称
samples <- c("anterior1", "posterior1")
#实际项目中需要将 LoadData() 替换为Load10X_Spatial()读取数据 。
stlist <- lapply(samples,function(sample_name) {
st_obj <- LoadData("stxBrain", type = sample_name)
st_obj$sample_id <- sample_name
SCTransform(st_obj, assay = "Spatial", verbose = FALSE)
})
names(stlist) <- samples
st <- merge(x = stlist[[1]],y = stlist[-1],
add.cell.ids = names(stlist))
DefaultAssay(st) <- "SCT"
3. 降维、整合、聚类分群
Harmony用于校正样本间批次效应，还有其他算法，且听下回分解。
# 每个样本单独计算高变基因，再取交集
features <- Reduce(intersect, lapply(stlist, VariableFeatures))
VariableFeatures(st) <- features
rm(stlist);invisible(gc())
# PCA
st <- RunPCA(st, verbose = FALSE)
#维度数量
dims <- 1:30
st <- RunHarmony(st,
group.by.vars = "sample_id",
dims.use = dims,
verbose = FALSE)
st <- FindNeighbors(st,
reduction = "harmony",
dims = dims,
verbose = FALSE)
st <- FindClusters(st, resolution = 0.5, verbose = FALSE)
st <- RunUMAP(st,
reduction = "harmony",
dims = dims,
verbose = FALSE)
UMAP图展示聚类结果和样本来源，SpatialDimPlot展示每个Cluster的空间分布。
p1 <- DimPlot(st, group.by = "sample_id") +
ggtitle("Samples")+
coord_fixed()
p2 <- DimPlot(st,group.by = "seurat_clusters",
label = TRUE,repel = TRUE) +
ggtitle("Clusters")+
coord_fixed()
p1 + p2
SpatialDimPlot(st,
group.by = "seurat_clusters",
label = TRUE,label.size = 3,
ncol = length(Images(st)))
查看感兴趣的基因的空间分布
DefaultAssay(st) <- "SCT"
SpatialFeaturePlot(st,features = c("Hpca", "Plp1"),
ncol = length(Images(st)),
alpha = c(0.1, 1))
4. scRNA参考数据
这是小鼠皮层细胞的高质量参考数据，来自 Allen 研究所。下载自 https://www.dropbox.com/s/cuowvm4vrf65pvq/allen_cortex.rds?dl=1
设置 ncells=3000 会将整个数据集归一化，但仅在 3000 个细胞上学习噪声模型。 这能显著加快 SCTransform 的速度，且性能无损失。是官网推荐的用法。
allen_reference <- readRDS("allen_cortex.rds")
dim(allen_reference)
## [1] 34617 14249
allen_reference <- SCTransform(allen_reference,
ncells = 3000,
verbose = FALSE) %>%
RunPCA(verbose = FALSE) %>%
RunUMAP(dims = 1:30)
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 | 血管及软脑膜细胞 |
5. 标签转移
寻找锚点，进行标签转移。这里的 st 已经完成了 SCTransform 和 PCA，因此可以直接用于寻找锚点和标签转移。
anchors <- FindTransferAnchors(reference = allen_reference,
query = st,
normalization.method = "SCT")
predictions.assay <- TransferData(anchorset = anchors,
refdata = allen_reference$subclass,
prediction.assay = TRUE,
weight.reduction = st[["harmony"]],
dims = dims)
GetAssayData(predictions.assay, layer = "data")[1:4, 1:4]
## anterior1_AAACAAGTATCTCCCA-1 anterior1_AAACACCAATAACTGC-1
## Vip 0 0
## Lamp5 0 0
## Sst 0 0
## Sncg 0 0
## anterior1_AAACAGAGCGACTCCT-1 anterior1_AAACAGCTTTCAGAAG-1
## Vip 0 0
## Lamp5 0 0
## Sst 0 0
## Sncg 0 0
此时我们得到的predictions.assay是每个spot是每一种细胞的预测分数。把他插入空转Seurat对象作为一个组成部分。
st[["predictions"]] <- predictions.assay
DefaultAssay(st) <- "predictions"
画图查看其中两种细胞的预测分数。
features <- c("Astro","L2/3 IT")
SpatialFeaturePlot(st,
features = features,
pt.size.factor = 1.6,
ncol = length(Images(st)),
alpha = c(0.1, 1))
以第一个样本为例，寻找(预测分数)有明显空间分布模式的细胞类型。注意，FindSpatiallyVariableFeatures是应该每个样本单独计算的。切换第二个样本时，把 image = samples[1]里面的1改为2即可。
st <- FindSpatiallyVariableFeatures(st,
assay = "predictions",
layer = "data",
features = setdiff(rownames(st), "max"),
image = samples[1],
selection.method = "moransi")
展示 Moran’s I 最高的 4 种细胞类型，也就是预测分数空间自相关性最强、空间分布模式最明显的细胞类型。（图上画的是细胞类型预测分数，不是基因表达量，也不是真实细胞数量或者丰度。）。
top.clusters <- st %>%
SpatiallyVariableFeatures(method = "moransi",
assay = "predictions") %>%
head(4)
top.clusters
## [1] "L6 CT" "Meis2" "L4" "Macrophage"
SpatialFeaturePlot(object = st, features = top.clusters, ncol = length(Images(st))) &
theme(plot.margin = unit(c(1, 1, 1, 1), "mm"),
legend.text = element_text(size = 8))
6. 保存结果
最终对象包含 PCA、Harmony embedding、联合聚类以及每个 spot 完整的细胞类型预测分数。
saveRDS(st, "st.rds")
转载自 CSDN-专业IT技术社区
原文链接：https://blog.csdn.net/weixin_42960896/article/details/163558649
