单细胞空间转录分析之Seurat-多样本整合（浅谈空间批次）

单细胞空间转录组和单细胞分析类似，不可避免的会遇到多样本的问题，这就需要使用多样本整合分析策略
Seurat提供了多张切片（slices）整合分析(Merge合并)

Seurat

单细胞空间转录分析之Seurat：https://www.jianshu.com/p/c9a601ced91f
单细胞空间转录分析之Seurat-多样本整合（浅谈空间批次）：https://www.jianshu.com/p/609b04096b79

从10X官网，我们可以下载获得小鼠大脑同一张切片前端和后端空转数据，包括anterior1，anterior2，posterior1，posterior2。https://support.10xgenomics.com/spatial-gene-expression/datasets/1.1.0/V1_Mouse_Brain_Sagittal_Anterior
https://support.10xgenomics.com/spatial-gene-expression/datasets/1.1.0/V1_Mouse_Brain_Sagittal_Anterior_Section_2
https://support.10xgenomics.com/spatial-gene-expression/datasets/1.1.0/V1_Mouse_Brain_Sagittal_Posterior
https://support.10xgenomics.com/spatial-gene-expression/datasets/1.1.0/V1_Mouse_Brain_Sagittal_Posterior_Section_2
导入包

library(Seurat)
library(ggplot2)
library(patchwork)
library(cowplot)
library(dplyr)

读取数据，并分别进行SCT标准化

dir = c('./Mouse/Brain_Section1_Sagittal_Anterior/Brain_anterior1/outs',
        './Mouse/Brain_Section2_Sagittal_Anterior/Brain_anterior2/outs',
        './Mouse/Brain_Section1_Sagittal_Posterior/Brain_posterior1/outs',
        './Mouse/Brain_Section2_Sagittal_Posterior/Brain_posterior2/outs')
names(dir) = c('anterior1', 'anterior2', 'posterior1', 'posterior2')
brain <- list()
for(i in 1:length(dir)){
brain[[i]] <-Seurat::Load10X_Spatial(data.dir = dir[i])
brain[[i]]@meta.data$orig.ident <-names(dir)[i]
}
for (i in 1:length(brain)) {
    brain[[i]] <- SCTransform(brain[[i]], assay = "Spatial", return.only.var.genes = FALSE, verbose = FALSE)} ##SCT标准化

merge合并

setwd("./Merge_Data") #设置输出路径
brain.merge <- merge(brain[[1]], y=c(brain[[2]], brain[[3]], brain[[4]]))
dim(brain.merge)
table(brain.merge@meta.data$orig.ident)
DefaultAssay(brain.merge) <- "SCT"
VariableFeatures(brain.merge) <- c(VariableFeatures(brain[[1]]), VariableFeatures(brain[[2]]), VariableFeatures(brain[[3]]), VariableFeatures(brain[[4]])) 

降维，聚类，可视化 这儿和单一样本处理方式相同

brain.merge <- RunPCA(brain.merge, verbose = FALSE)
brain.merge <- FindNeighbors(brain.merge, dims = 1:30)
brain.merge <- FindClusters(brain.merge,resolution = 0.8, verbose = FALSE)
brain.merge <- RunUMAP(brain.merge, dims = 1:30)
pdf("merge_umap.pdf",width=10,height=5)
DimPlot(brain.merge, reduction = "umap", group.by = c("ident", "orig.ident"))
dev.off()
pdf("merge_umap_slice.pdf",width=10,height=5)
SpatialDimPlot(brain.merge) 
dev.off()
pdf("merge_umap_gene.pdf",width=10,height=5)
SpatialFeaturePlot(brain.merge, features = c("Hpca", "Plp1"))
dev.off()

merge_umap
image.png
merge_umap_gene
识别簇marker genes

saveRDS(brain.merge,"merge.rds") #保存数据
brain.markers <- FindAllMarkers(brain.merge, only.pos = FALSE, min.pct = 0.25, logfc.threshold = 0.25,test.use = "wilcox")
write.table(brain.markers,"marker.txt",row.names=TRUE,col.names=TRUE,sep="\t")
topgene<-brain.markers %>% group_by(cluster) %>% top_n(n = 3, wt = avg_logFC)
pdf("merge_heatmap.pdf",width=10,height=8)
DoHeatmap(brain.merge, features = topgene$gene,size = 2) + NoLegend()
dev.off()

image.png

我们可以发现重复切片之间分类几乎一样，但是anterior和posterior之间有很大差异，尤其中间切断位置不能连接起来，观察marker genes表达，进一步说明merge直接整合，分类不是很准确，样本间还是有批次存在的，这儿我们可以使用单细胞去批次方法去除空间批次。

这儿，我们可以尝试使用Seurat FindIntegrationAnchors函数来整合表达谱，去除批次。
读取数据

dir = c('./Mouse/Brain_Section1_Sagittal_Anterior/Brain_anterior1/outs',
        './Mouse/Brain_Section1_Sagittal_Posterior/Brain_posterior1/outs')
names(dir) = c('anterior1', 'posterior1')
brain <- list()
for(i in 1:length(dir)){
brain[[i]] <-Seurat::Load10X_Spatial(data.dir = dir[i])
brain[[i]]@meta.data$orig.ident <-names(dir)[i]
}

标准化NormalizeData

for (i in 1:length(brain)) {
    brain[[i]] <- NormalizeData(brain[[i]]) #未使用SCT标准化
    brain[[i]] <- FindVariableFeatures(brain[[i]], selection.method = "vst", nfeatures = 2000)}
##
features <- SelectIntegrationFeatures(object.list = brain)
brain.anchors <- FindIntegrationAnchors(object.list = brain, anchor.features = features)
brain_in <- IntegrateData(anchorset = brain.anchors)
dim(brain_in)
DefaultAssay(brain_in) <- "integrated"
brain_in <- ScaleData(brain_in, verbose = FALSE)
brain_in <- RunPCA(brain_in, npcs = 30, verbose = FALSE)
brain_in <- RunUMAP(brain_in, reduction = "pca", dims = 1:30)
brain_in <- FindNeighbors(brain_in, reduction = "pca", dims = 1:30)
brain_in <- FindClusters(brain_in, resolution = 0.8,verbose = FALSE)
pdf("merge_umap1.pdf",width=10,height=5)
DimPlot(brain_in, reduction = "umap", group.by = c("ident", "orig.ident"))
dev.off()
pdf("merge_umap_slice1.pdf",width=10,height=5)
SpatialDimPlot(brain_in)
dev.off()

merge_umap1.
merge_umap_slice1
图片上看，两张切片连接处分类相比未去批次前更为合理。
SCT标准化

for (i in 1:length(brain)) {
    brain[[i]] <- SCTransform(brain[[i]], assay = "Spatial", return.only.var.genes = FALSE, verbose = FALSE)} #使用SCT标准化
features <- SelectIntegrationFeatures(object.list = brain, nfeatures = 3000)
brain <- PrepSCTIntegration(object.list = brain, anchor.features = features)
brain.anchors <- FindIntegrationAnchors(object.list = brain, normalization.method = "SCT", anchor.features = features)
brian.sct <- IntegrateData(anchorset = brain.anchors, normalization.method = "SCT")
dir.create("Integrate_SCT_data")
setwd("Integrate_SCT_data")
dim(brian.sct)
brian.sct <- RunPCA(brian.sct, npcs = 30, verbose = FALSE)
brian.sct <- RunUMAP(brian.sct, reduction = "pca", dims = 1:30)
brian.sct <- FindNeighbors(brian.sct, reduction = "pca", dims = 1:30)
brian.sct <- FindClusters(brian.sct, resolution = 0.8,verbose = FALSE)
pdf("merge_umap2.pdf",width=10,height=5)
DimPlot(brian.sct, reduction = "umap", group.by = c("ident", "orig.ident"))
dev.off()
pdf("merge_umap_slice2.pdf",width=10,height=5)
SpatialDimPlot(brian.sct)
dev.off()

merge_umap2
merge_umap_slice2

相比上面，感觉spots会混乱一些，一些类别中含有更多错误分类的的spots。