此教程是使用Seurat分析 Slide-seq v2數(shù)據(jù)；大部分流程和步驟還是和10x Visium分析流程一樣 。

# 數(shù)據(jù)集

InstallData("ssHippo")

slide.seq <- LoadData("ssHippo")

# 數(shù)據(jù)預(yù)處理

這一步與10x Visium數(shù)據(jù)處理差不多；但是當(dāng)前的測試數(shù)據(jù)含有許多UMI很低的值，但在本教程中也是保留所有檢測到的數(shù)據(jù)進(jìn)行下游分析。

plot1 <- VlnPlot(slide.seq, features = "nCount_Spatial", pt.size = 0, log = TRUE) + NoLegend()
slide.seq$log_nCount_Spatial <- log(slide.seq$nCount_Spatial)
plot2 <- SpatialFeaturePlot(slide.seq, features = "log_nCount_Spatial") + theme(legend.position = "right")
wrap_plots(plot1, plot2)

image.png

然后使用sctransform對數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化，并進(jìn)行數(shù)據(jù)降維和聚類。

slide.seq <- SCTransform(slide.seq, assay = "Spatial", ncells = 3000, verbose = FALSE)
slide.seq <- RunPCA(slide.seq)
slide.seq <- RunUMAP(slide.seq, dims = 1:30)
slide.seq <- FindNeighbors(slide.seq, dims = 1:30)
slide.seq <- FindClusters(slide.seq, resolution = 0.3, verbose = FALSE)

基于UMAP 可視化數(shù)據(jù)結(jié)果（DimPlot()）或給予bead對應(yīng)位置進(jìn)行空間可視化（SpatialDimPlot()）

plot1 <- DimPlot(slide.seq, reduction = "umap", label = TRUE)
plot2 <- SpatialDimPlot(slide.seq, stroke = 0)
plot1 + plot2

image.png

SpatialDimPlot(slide.seq, cells.highlight = CellsByIdentities(object = slide.seq, idents = c(1,
    6, 13)), facet.highlight = TRUE)

image.png

# 整合scRNA-seq參考數(shù)據(jù)

為了便于對Slide seq數(shù)據(jù)集進(jìn)行細(xì)胞類型注釋，利用現(xiàn)有的小鼠ScRNAseq海馬數(shù)據(jù)集( Saunders, Macosko, et al. 2018)。

處理后的數(shù)據(jù)，保存為Seurat對象，下載地址： data

原始矩陣數(shù)據(jù)下載： DropViz website

ref <- readRDS("./mouse_hippocampus_reference.rds")

paper中細(xì)胞類型注釋保存在Seurat 的cell metadata中；注釋包含了不同層級的注釋，ref$class是初始注釋，ref$subcluster是對初始注釋進(jìn)行了的下一級注釋；

首先運(yùn)行Seurat標(biāo)簽轉(zhuǎn)移方法來預(yù)測每個磁柱上的主要細(xì)胞類型。

anchors <- FindTransferAnchors(reference = ref, query = slide.seq, normalization.method = "SCT",
    npcs = 50)
predictions.assay <- TransferData(anchorset = anchors, refdata = ref$celltype, prediction.assay = TRUE,
    weight.reduction = slide.seq[["pca"]], dims = 1:50)
slide.seq[["predictions"]] <- predictions.assay

然后可視化預(yù)測的分?jǐn)?shù)：

DefaultAssay(slide.seq) <- "predictions"
SpatialFeaturePlot(slide.seq, features = c("Dentate Principal cells", "CA3 Principal cells", "Entorhinal cortex",
    "Endothelial tip", "Ependymal", "Oligodendrocyte"), alpha = c(0.1, 1))

slide.seq$predicted.id <- GetTransferPredictions(slide.seq)
Idents(slide.seq) <- "predicted.id"
SpatialDimPlot(slide.seq, cells.highlight = CellsByIdentities(object = slide.seq, idents = c("CA3 Principal cells",
    "Dentate Principal cells", "Endothelial tip")), facet.highlight = TRUE)

image.png

# 鑒定空間差異表達(dá)基因

一般使用兩種方法：

• 統(tǒng)計分析預(yù)注釋的切片區(qū)域鑒定差異表達(dá)基因；
• 統(tǒng)計分析基因表達(dá)對空間位置的依賴性

對于第二種方法，FindSpatiallyVariableFeatures(method = 'moransi')便可以實(shí)現(xiàn)，函數(shù)實(shí)現(xiàn)了一種叫Moran’s I 的算法；Moran’s I 可以對一個基因表達(dá)對空間位置的依賴性進(jìn)行打分（類似于相關(guān)性系數(shù)）。這就允許我們基于基因的空間表達(dá)變化進(jìn)行打分。這就允許我們依賴于基因表達(dá)值對基因進(jìn)行快速的排序。為了快速實(shí)現(xiàn)這個計算過程，seurat采用了一個基本區(qū)間(bin)的策略，基于選擇一個矩陣方框的Slide-seq puck，然后平均其表達(dá)和定位。在x和y軸上bin的數(shù)量分別使用x.cuts和y.cuts確定。另外，可以安裝Rfast2包install.packages('Rfast2'), 這個包可以減少這個過程運(yùn)行時間。

DefaultAssay(slide.seq) <- "SCT"
slide.seq <- FindSpatiallyVariableFeatures(slide.seq, assay = "SCT", slot = "scale.data", features = VariableFeatures(slide.seq)[1:1000],
    selection.method = "moransi", x.cuts = 100, y.cuts = 100)

展示Moran’s I鑒定的top6個基因：

SpatialFeaturePlot(slide.seq, features = head(SpatiallyVariableFeatures(slide.seq, selection.method = "moransi"),
    6), ncol = 3, alpha = c(0.1, 1), max.cutoff = "q95")

image.png

# 原文
Analysis, visualization, and integration of spatial datasets with Seurat

# 系列文章
Seurat分析10x Visium空間轉(zhuǎn)錄組數(shù)據(jù)
單細(xì)胞測序分析: Seurat V3聯(lián)合harmony進(jìn)行單細(xì)胞數(shù)據(jù)整合分析
單細(xì)胞測序分析: Seurat 使用教程