单细胞空间转录分析之Seurat

将空间位置信息和转录组分析相结合，对于癌症、免疫、肿瘤免疫相互作用，组织微环境，神经和发育等领域，有着令人期待的应用前景。

而作为单细胞转录组必备R包Seurat,与时俱进，可用来分析单细胞空间转录组数据。https://satijalab.org/seurat/articles/spatial_vignette.html

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

和分析单细胞转录组数据一样，单细胞空间转录组主要包括了：质控（QC），标准化（Normalization），降维聚类（Dimensional reduction and clustering），Cluster marker genes， Spatially variable genes

这儿使用10X官网提供的单细胞转录组数据集：https://support.10xgenomics.com/spatial-gene-expression/datasets/1.1.0/V1_Mouse_Brain_Sagittal_Anterior， 新鲜的冷冻小鼠脑组织， 前牙矢状切面，可以参考前面讲述的ABA大脑图谱：https://www.jianshu.com/p/5d087fffeb35

导入相关包

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

读取数据

brain <-Seurat::Load10X_Spatial(data.dir = "/home/wucheng/data_set/Spatial/Mouse/Brain_Section1_Sagittal_Anterior/Brain_anterior1/outs")
dir.create("/home/wucheng/jianshu/seurat/spatial")  #创建输出路径
setwd("/home/wucheng/jianshu/seurat/spatial") #提供文件输出路径
brain  #查看包含的spots数和基因数，S4格式
An object of class Seurat 
32285 features across 2695 samples within 1 assay 
Active assay: Spatial (32285 features, 0 variable features)

质控 QC 查看总counts和表达的genes在spots中的分布：

plot1 <- VlnPlot(brain, features = "nCount_Spatial", pt.size = 0.1) + NoLegend()
plot2 <- VlnPlot(brain, features = "nFeature_Spatial", pt.size = 0.1) + NoLegend()
pdf("QC.pdf",width=10,height=5) 
print(wrap_plots(plot1, plot2))
dev.off()

QC
标准化 Seurat提出 LogNormalize函数可能会有问题，因为它会强制每个数据点在标准化之后具有相同的底层“大小”,替代方法，推荐使用sctransform

brain <- SCTransform(brain, assay = "Spatial", return.only.var.genes = FALSE, verbose = FALSE) #包含NormalizeData, FindVariableFeatures, ScaleData workflow 一步解决预处理

降维聚类可视化 Dimensionality reduction, clustering, and visualization，和单细胞转录组分析一致

brain <- RunPCA(brain, assay = "SCT", verbose = FALSE)
brain<- FindNeighbors(brain, reduction = "pca", dims = 1:30)
brain<- FindClusters(brain, resolution = 0.8, verbose = FALSE)  ##resolution分辨率可改变，默认0.8
brain <- RunUMAP(brain, reduction = "pca", dims = 1:30)
p1 <- FeaturePlot(brain, features = "nCount_Spatial")  + theme(legend.position = "right")
p2 <- DimPlot(brain, reduction = "umap", label = TRUE)
pdf("umap_cluster.pdf",width=12,height=5)
print(plot_grid(p1, p2))
dev.off()
p1 <- SpatialFeaturePlot(brain, features = "nCount_Spatial") + theme(legend.position = "right")
p2 <- SpatialDimPlot(brain, label = TRUE, label.size = 3)
pdf("spatial_cluster.pdf",width=12,height=5)
print(plot_grid(p1, p2))
dev.off()

umap_Cluster
spatial_Cluster
显示每簇位置

clu<-length(unique(brain@ active.ident))
pdf("spatial_cluster_sub.pdf",width=6,height=12)
print(SpatialDimPlot(brain, cells.highlight = CellsByIdentities(object = brain, idents = c(0:(clu-1))), facet.highlight = TRUE, ncol = 3))
dev.off()

spatial_cluster_sub
关键基因的表达可视化 #Hpca is a strong hippocampus marker and Ttr is a marker of the choroid plexus

spatial_gene_exp
spatial_gene_exp1

Seurat提供了两种工作流程来识别与组织内空间位置相关的分子特征。第一种是基于组织内预先标注的解剖区域执行差异表达，这可以从无监督的聚类或先验知识中确定。在这种情况下，此策略将起作用，因为上面的群集显示出明显的空间限制。
每一簇marker genes

brain.markers <- FindAllMarkers(brain, only.pos = FALSE, min.pct = 0.25, logfc.threshold = 0.25,test.use = "wilcox")
head(brain.markers)
                 p_val avg_logFC pct.1 pct.2     p_val_adj cluster     gene
Ido1     3.922047e-286  0.881875 0.757 0.058 7.165579e-282       0     Ido1
Drd2     2.402721e-280  1.551381 0.965 0.206 4.389771e-276       0     Drd2
Lrrc10b  1.745108e-262  1.481596 0.965 0.240 3.188313e-258       0  Lrrc10b
Adora2a  2.005073e-258  1.588071 0.981 0.274 3.663268e-254       0  Adora2a
Gpr6     2.482538e-254  1.167760 0.887 0.148 4.535597e-250       0     Gpr6
Syndig1l 3.530456e-250  1.517816 0.986 0.321 6.450144e-246       0 Syndig1l

write.table(brain.markers,"marker.txt",row.names=TRUE,col.names=TRUE,sep="\t")
de_markers <- FindMarkers(brain, ident.1 = 5, ident.2 = 6) #或者识别任意两簇的差异
head(de_markers)
               p_val avg_logFC pct.1 pct.2    p_val_adj
Calb2   3.000452e-68  2.326368 1.000 0.558 5.481826e-64
Camk2n1 4.835827e-68 -1.681583 1.000 1.000 8.835055e-64
Nrgn    4.983341e-68 -2.207203 0.981 1.000 9.104565e-64
Stx1a   8.443715e-68 -1.548973 0.797 1.000 1.542667e-63
Nptxr   2.764752e-67 -1.340551 0.942 1.000 5.051201e-63
Hpca    8.214733e-67 -1.597642 0.816 1.000 1.500832e-62
pdf("spatial_gene_exp_C5VSC6.pdf",width=15,height=5)
SpatialFeaturePlot(object = brain, features = rownames(de_markers)[1:3], alpha = c(0.1, 1), ncol = 3)
dev.off()

spatial_gene_exp_C5VSC6
另一种方法是在没有预先注释的情况下搜索表现出空间图案的特征，使用FindSpatiallyVariables,默认方法（method ='markvariogram）。
空间特异性genes

brain <- FindSpatiallyVariableFeatures(brain, assay = "SCT", features = VariableFeatures(brain)[1:1000],
    selection.method = "markvariogram")
top.features <- head(SpatiallyVariableFeatures(brain, selection.method = "markvariogram"), 20)
write.table(SpatiallyVariableFeatures(brain, selection.method = "markvariogram"),"SpatiallyVariableFeatures.txt")
pdf("Spatial_Var_feature.pdf",width=12,height=15)
print(SpatialFeaturePlot(brain, features = top.features, ncol = 4, alpha = c(0.1, 1)))
dev.off()
}

Spatial_Var_feature_top20

**保存meta文件和数据

write.table(brain@meta.data,"meta.data.txt",row.names=TRUE,col.names=TRUE,sep="\t")
saveRDS(brain,"data.rds")