將空間位置信息和轉(zhuǎn)錄組分析相結(jié)合堕绩，對(duì)于癌癥、免疫邑时、腫瘤免疫相互作用奴紧，組織微環(huán)境，神經(jīng)和發(fā)育等領(lǐng)域晶丘，有著令人期待的應(yīng)用前景黍氮。

而作為單細(xì)胞轉(zhuǎn)錄組必備R包Seurat,與時(shí)俱進(jìn)，可用來分析單細(xì)胞空間轉(zhuǎn)錄組數(shù)據(jù)浅浮。https://satijalab.org/seurat/articles/spatial_vignette.html

單細(xì)胞空間轉(zhuǎn)錄分析之Seurat：http://www.reibang.com/p/c9a601ced91f
單細(xì)胞空間轉(zhuǎn)錄分析之Seurat-多樣本整合（淺談空間批次）：http://www.reibang.com/p/609b04096b79

和分析單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)一樣沫浆，單細(xì)胞空間轉(zhuǎn)錄組主要包括了：質(zhì)控（QC），標(biāo)準(zhǔn)化（Normalization）滚秩，降維聚類（Dimensional reduction and clustering）件缸，Cluster marker genes， Spatially variable genes

這兒使用10X官網(wǎng)提供的單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)集：https://support.10xgenomics.com/spatial-gene-expression/datasets/1.1.0/V1_Mouse_Brain_Sagittal_Anterior叔遂， 新鮮的冷凍小鼠腦組織他炊， 前牙矢狀切面，可以參考前面講述的ABA大腦圖譜：http://www.reibang.com/p/5d087fffeb35

導(dǎo)入相關(guān)包

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

讀取數(shù)據(jù)

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")  #創(chuàng)建輸出路徑
setwd("/home/wucheng/jianshu/seurat/spatial") #提供文件輸出路徑
brain  #查看包含的spots數(shù)和基因數(shù)已艰，S4格式
An object of class Seurat 
32285 features across 2695 samples within 1 assay 
Active assay: Spatial (32285 features, 0 variable features)

質(zhì)控 QC 查看總counts和表達(dá)的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

brain <- SCTransform(brain, assay = "Spatial", return.only.var.genes = FALSE, verbose = FALSE) #包含NormalizeData, FindVariableFeatures, ScaleData workflow 一步解決預(yù)處理

降維聚類可視化 Dimensionality reduction, clustering, and visualization哩掺，和單細(xì)胞轉(zhuǎn)錄組分析一致

brain <- RunPCA(brain, assay = "SCT", verbose = FALSE)
brain<- FindNeighbors(brain, reduction = "pca", dims = 1:30)
brain<- FindClusters(brain, resolution = 0.8, verbose = FALSE)  ##resolution分辨率可改變凿叠，默認(rèn)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

spatial_gene_exp

spatial_gene_exp1

Seurat提供了兩種工作流程來識(shí)別與組織內(nèi)空間位置相關(guān)的分子特征。第一種是基于組織內(nèi)預(yù)先標(biāo)注的解剖區(qū)域執(zhí)行差異表達(dá)，這可以從無監(jiān)督的聚類或先驗(yàn)知識(shí)中確定盒件。在這種情況下蹬碧，此策略將起作用，因?yàn)樯厦娴娜杭@示出明顯的空間限制炒刁。
每一簇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) #或者識(shí)別任意兩簇的差異
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

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文件和數(shù)據(jù)

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