第十七計 拋磚引玉
以自己的粗淺的意見引出別人高明的見解惊科。

概述

本教程演示了如何使用Seurat（> = 3.2）分析空間分辨的RNA-seq數(shù)據(jù)。盡管分析流程類似于單細胞RNA序列分析的Seurat工作流程蒲跨，但我們引入了更新的交互作用和可視化工具译断，特別著重于空間和分子信息的整合授翻。本教程將涵蓋以下任務(wù)或悲，我們認為這些任務(wù)在許多空間分析中都是常見的：

• 正常化
• 降維和聚類
• 檢測空間可變特征
• 互動可視化
• 與單細胞RNA-seq數(shù)據(jù)整合
• 處理多個切片

對于我們的第一個小插圖堪唐，我們分析了使用10x Genomics的Visium技術(shù)生成的數(shù)據(jù)集巡语。我們將擴展Seurat以便在不久的將來使用其他數(shù)據(jù)類型，包括SLIDE-Seq淮菠，STARmap和MERFISH男公。

首先，我們加載Seurat和此小插圖所需的其他軟件包合陵。

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

10倍的Visium

數(shù)據(jù)集

在這里枢赔，我們將使用使用Visium v1化學方法生成的最近發(fā)布的矢狀小鼠大腦切片的數(shù)據(jù)集澄阳。有兩個連續(xù)的前節(jié)和兩個（匹配的）連續(xù)后節(jié)。

您可以在此處下載數(shù)據(jù)踏拜，然后使用該Load10X_Spatial()功能將其加載到Seurat中碎赢。這將讀取spaceranger管道的輸出，并返回一個Seurat對象速梗，該對象包含點級表達數(shù)據(jù)以及組織切片的關(guān)聯(lián)圖像肮塞。您還可以使用我們的SeuratData包來輕松訪問數(shù)據(jù)，如下所示姻锁。安裝數(shù)據(jù)集后枕赵，您可以鍵入?stxBrain以了解更多信息。

InstallData("stxBrain")

brain <- LoadData("stxBrain", type = "anterior1")

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

我們通過基因表達數(shù)據(jù)當場執(zhí)行的初始預處理步驟與典型的scRNA-seq實驗相似位隶。我們首先需要對數(shù)據(jù)進行歸一化拷窜，以解決各個數(shù)據(jù)點之間的測序深度差異。我們注意到钓试，對于空間數(shù)據(jù)集装黑，分子計數(shù)/斑點的差異可能很大，尤其是在組織中細胞密度存在差異的情況下弓熏。我們在這里看到大量的異質(zhì)性恋谭，這需要有效的規(guī)范化。

plot1 <- VlnPlot(brain, features = "nCount_Spatial", pt.size = 0.1) + NoLegend()
plot2 <- SpatialFeaturePlot(brain, features = "nCount_Spatial") + theme(legend.position = "right")
wrap_plots(plot1, plot2)

image.png

這些圖表明挽鞠，斑點上分子計數(shù)的變化不僅是技術(shù)上的問題疚颊，而且還取決于組織的解剖結(jié)構(gòu)。例如信认，神經(jīng)元（例如皮質(zhì)白質(zhì)）耗竭的組織區(qū)域可再現(xiàn)地顯示出較低的分子數(shù)材义。結(jié)果，在標準化后LogNormalize()`強制每個數(shù)據(jù)點具有相同的底層“大小”的標準方法（例如函數(shù)）可能會出現(xiàn)問題嫁赏。

作為替代方案其掂，我們建議使用sctransform（Hafemeister和Satija，Genome Biology 2019）潦蝇，它構(gòu)建了基因表達的正則化負二項式模型款熬，以便在保留生物學差異的同時考慮技術(shù)偽像朵逝。有關(guān)sctransform更多詳細信息庇忌，請參見文章在這里和修拉的小插曲在這里。sctransform可以對數(shù)據(jù)進行歸一化王暗，檢測高變異特征并將數(shù)據(jù)存儲在SCT測定中则酝。

brain <- SCTransform(brain, assay = "Spatial", verbose = FALSE)

基因表達可視化

在Seurat中殉簸，我們具有探索空間數(shù)據(jù)固有的視覺本質(zhì)并與之交互的功能。SpatialFeaturePlot()Seurat中的功能FeaturePlot()`可以擴展，并且可以在組織組織學之上疊加分子數(shù)據(jù)般卑。例如武鲁，在此小鼠大腦數(shù)據(jù)集中，基因Hpca是強海馬標志物蝠检，而Ttr是脈絡(luò)叢的標志物洞坑。

SpatialFeaturePlot(brain, features = c("Hpca", "Ttr"))

image.png

Seurat中的默認參數(shù)強調(diào)分子數(shù)據(jù)的可視化。但是蝇率，您還可以通過更改以下參數(shù)來調(diào)整斑點的大谐僭印（及其透明度），以改善組織學圖像的可視化：

• pt.size.factor-這將縮放斑點的大小本慕。默認值為1.6
• alpha-最小和最大透明度排拷。默認值為c（1，1）锅尘。
• 嘗試設(shè)置為alphac（0.1监氢，1），以降低具有較低表達式的點的透明度

p1 <- SpatialFeaturePlot(brain, features = "Ttr", pt.size.factor = 1)
p2 <- SpatialFeaturePlot(brain, features = "Ttr", alpha = c(0.1, 1))
p1 + p2

image

降維藤违，聚類和可視化

然后浪腐，我們可以使用與scRNA-seq分析相同的工作流程，對RNA表達數(shù)據(jù)進行降維和聚類顿乒。

brain <- RunPCA(brain, assay = "SCT", verbose = FALSE)
brain <- FindNeighbors(brain, reduction = "pca", dims = 1:30)
brain <- FindClusters(brain, verbose = FALSE)
brain <- RunUMAP(brain, reduction = "pca", dims = 1:30)

然后议街，我們可以在UMAP空間帶有DimPlot()）中可視化聚類的結(jié)果，或者使用覆蓋可視化圖像上的聚類結(jié)果SpatialDimPlot()

p1 <- DimPlot(brain, reduction = "umap", label = TRUE)
p2 <- SpatialDimPlot(brain, label = TRUE, label.size = 3)
p1 + p2

image

由于存在許多種顏色璧榄，因此可視化哪個體素屬于哪個群集可能具有挑戰(zhàn)性特漩。我們有一些策略可以幫助您解決這個問題。設(shè)置label參數(shù)會在每個群集的中位數(shù)處放置一個彩色框（請參見上圖）骨杂。

您還可以使用cells.highlight參數(shù)在上劃分特定的關(guān)注單元格SpatialDimPlot()`涂身。如下所示，這對于區(qū)分單個群集的空間定位非常有用搓蚪。

SpatialDimPlot(brain, cells.highlight = CellsByIdentities(object = brain, idents = c(2, 1, 4, 3, 
    5, 8)), facet.highlight = TRUE, ncol = 3)

image

交互式繪圖

我們還內(nèi)置了許多交互式繪圖功能蛤售。無論SpatialDimPlot()和SpatialFeaturePlot()現(xiàn)在有一個interactive參數(shù)，當設(shè)置為TRUE妒潭，將打開Rstudio觀眾面板與互動閃亮的情節(jié)悴能。下面的示例演示了一種交互式SpatialDimPlot()的方法，您可以在其上懸停并查看單元名稱和當前標識類（類似于先前的do.hover`行為）杜耙。

SpatialDimPlot(brain, interactive = TRUE)

對于SpatialFeaturePlot()搜骡，將“交互性”設(shè)置為TRUE會彈出一個交互窗格拂盯，您可以在其中調(diào)整點的透明度佑女，點的大小以及Assay`要繪制的和特征。瀏覽數(shù)據(jù)后，選擇完成按鈕將返回最后一個活動圖作為ggplot對象团驱。

SpatialFeaturePlot(brain, features = "Ttr", interactive = TRUE)

該LinkedDimPlot()`功能將UMAP表示鏈接到組織圖像表示摸吠，并允許交互式選擇。例如嚎花，您可以在UMAP圖中選擇一個區(qū)域寸痢，并且圖像表示中的相應斑點將突出顯示。

LinkedDimPlot(brain)

識別空間可變特征

Seurat提供了兩種工作流程來識別與組織內(nèi)空間位置相關(guān)的分子特征紊选。第一種是基于組織內(nèi)預先標注的解剖區(qū)域執(zhí)行差異表達啼止，這可以從無監(jiān)督的聚類或先驗知識中確定。在這種情況下兵罢，此策略將起作用献烦，因為上面的群集顯示出明顯的空間限制。

de_markers <- FindMarkers(brain, ident.1 = 5, ident.2 = 6)
SpatialFeaturePlot(object = brain, features = rownames(de_markers)[1:3], alpha = c(0.1, 1), ncol = 3)

image.png

在中實施的另一種方法FindSpatiallyVariables()是在沒有預先注釋的情況下搜索表現(xiàn)出空間圖案的特征卖词。默認方法（method = 'markvariogram）受Trendsceek啟發(fā)巩那，該工具將空間轉(zhuǎn)錄組學數(shù)據(jù)建模為標記點過程，并計算“變異函數(shù)”此蜈，該變異函數(shù)可識別表達水平取決于其空間位置的基因即横。更具體地說，此過程計算gamma（r）值裆赵，該值測量相距某個“ r”距離的兩個點之間的依賴性东囚。默認情況下，我們在這些分析中使用r值“ 5”战授，并且僅計算可變基因的這些值（其中變異獨立于空間位置而計算）以節(jié)省時間舔庶。

我們注意到，文獻中有多種方法可以完成此任務(wù)陈醒，包括SpatialDE和Splotch惕橙。我們鼓勵感興趣的用戶探索這些方法，并希望在不久的將來為其提供支持钉跷。

brain <- FindSpatiallyVariableFeatures(brain, assay = "SCT", features = VariableFeatures(brain)[1:1000], 
    selection.method = "markvariogram")

現(xiàn)在弥鹦，我們可視化此度量確定的前6個特征的表達。

top.features <- head(SpatiallyVariableFeatures(brain, selection.method = "markvariogram"), 6)
SpatialFeaturePlot(brain, features = top.features, ncol = 3, alpha = c(0.1, 1))

image.png

細分解剖區(qū)域

與單單元格對象一樣爷辙，您可以對對象進行子集化以集中處理數(shù)據(jù)的子集彬坏。在這里，我們大約將額葉皮層子集化膝晾。該過程還有助于在下一節(jié)中將這些數(shù)據(jù)與皮質(zhì)scRNA-seq數(shù)據(jù)集進行整合栓始。首先，我們獲取群集的子集血当，然后根據(jù)精確位置進行進一步細分幻赚。子集化后禀忆，我們可以在完整圖像或裁剪后的圖像上可視化皮質(zhì)細胞。

cortex <- subset(brain, idents = c(1, 2, 3, 4, 6, 7))
# now remove additional cells, use SpatialDimPlots to visualize what to remove
# SpatialDimPlot(cortex,cells.highlight = WhichCells(cortex, expression = image_imagerow > 400 |
# image_imagecol < 150))
cortex <- subset(cortex, anterior1_imagerow > 400 | anterior1_imagecol < 150, invert = TRUE)
cortex <- subset(cortex, anterior1_imagerow > 275 & anterior1_imagecol > 370, invert = TRUE)
cortex <- subset(cortex, anterior1_imagerow > 250 & anterior1_imagecol > 440, invert = TRUE)

p1 <- SpatialDimPlot(cortex, crop = TRUE, label = TRUE)
p2 <- SpatialDimPlot(cortex, crop = FALSE, label = TRUE, pt.size.factor = 1, label.size = 3)
p1 + p2

image

與單細胞數(shù)據(jù)集成

在?50um時落恼，來自于病毒測定的斑點將涵蓋多個細胞的表達譜箩退。對于可獲得scRNA-seq數(shù)據(jù)的系統(tǒng)列表不斷增加，用戶可能有興趣對每個空間體素進行“反卷積”以預測細胞類型的潛在組成佳谦。在準備此小插圖時戴涝，我們使用了參考scRNA-seq數(shù)據(jù)集測試了多種decovonlution和整合方法使用SMART-Seq2協(xié)議生成的來自Allen Institute的約14,000個成年小鼠皮質(zhì)細胞分類學。我們一直發(fā)現(xiàn)使用積分方法（而不是反卷積方法）具有更好的性能钻蔑，這可能是由于表征空間和單細胞數(shù)據(jù)集的噪聲模型存在很大差異啥刻，并且積分方法專門設(shè)計為對這些差異具有魯棒性。因此咪笑，我們應用了Seurat v3中引入的基于“錨”的集成工作流郑什，該工作流使注釋能夠從引用到查詢集的概率傳輸。因此蒲肋，我們利用sctransform歸一化方法遵循此處介紹的標簽傳輸工作流程蘑拯，但期望開發(fā)出新方法來完成此任務(wù)。

我們首先加載數(shù)據(jù)（可在此處下載）兜粘，對scRNA-seq參考進行預處理申窘，然后執(zhí)行標簽轉(zhuǎn)移。該過程為每個點輸出每個scRNA-seq派生類的概率分類孔轴。我們將這些預測添加為Seurat對象中的一項新分析剃法。

allen_reference <- readRDS("../data/allen_cortex.rds")

# note that setting ncells=3000 normalizes the full dataset but learns noise models on 3k cells
# this speeds up SCTransform dramatically with no loss in performance
library(dplyr)
allen_reference <- SCTransform(allen_reference, ncells = 3000, verbose = FALSE) %>% RunPCA(verbose = FALSE) %>% 
    RunUMAP(dims = 1:30)

# After subsetting, we renormalize cortex
cortex <- SCTransform(cortex, assay = "Spatial", verbose = FALSE) %>% RunPCA(verbose = FALSE)
# the annotation is stored in the 'subclass' column of object metadata
DimPlot(allen_reference, group.by = "subclass", label = TRUE)

image

anchors <- FindTransferAnchors(reference = allen_reference, query = cortex, normalization.method = "SCT")
predictions.assay <- TransferData(anchorset = anchors, refdata = allen_reference$subclass, prediction.assay = TRUE, 
    weight.reduction = cortex[["pca"]], dims = 1:30)
cortex[["predictions"]] <- predictions.assay

現(xiàn)在，我們獲得了每個班級每個地點的預測分數(shù)路鹰。在額葉皮層區(qū)域中特別令人感興趣的是層狀興奮性神經(jīng)元贷洲。在這里，我們可以區(qū)分這些神經(jīng)元亞型的不同順序?qū)咏纾?/p>

DefaultAssay(cortex) <- "predictions"
SpatialFeaturePlot(cortex, features = c("L2/3 IT", "L4"), pt.size.factor = 1.6, ncol = 2, crop = TRUE)

image

基于這些預測分數(shù)优构，我們還可以預測位置受空間限制的單元格類型。我們使用基于標記點過程的相同方法來定義空間可變特征雁竞，但將細胞類型預測得分用作“標記”而不是基因表達钦椭。

cortex <- FindSpatiallyVariableFeatures(cortex, assay = "predictions", selection.method = "markvariogram", 
    features = rownames(cortex), r.metric = 5, slot = "data")
top.clusters <- head(SpatiallyVariableFeatures(cortex), 4)
SpatialPlot(object = cortex, features = top.clusters, ncol = 2)

image.png

最后，我們證明我們的整合程序能夠恢復神經(jīng)元和非神經(jīng)元子集（包括層狀興奮性碑诉，第1層星形膠質(zhì)細胞和皮質(zhì)灰質(zhì)）的已知空間定位模式彪腔。

SpatialFeaturePlot(cortex, features = c("Astro", "L2/3 IT", "L4", "L5 PT", "L5 IT", "L6 CT", "L6 IT", 
    "L6b", "Oligo"), pt.size.factor = 1, ncol = 2, crop = FALSE, alpha = c(0.1, 1))

image.png

在Seurat中處理多個切片

小鼠大腦的該數(shù)據(jù)集包含與大腦另一半相對應的另一個切片。在這里进栽，我們將其讀入并執(zhí)行相同的初始歸一化德挣。

brain2 <- LoadData("stxBrain", type = "posterior1")
brain2 <- SCTransform(brain2, assay = "Spatial", verbose = FALSE)

為了在同一個Seurat對象中使用多個切片，我們提供了該merge功能快毛。

brain.merge <- merge(brain, brain2)

然后格嗅，這使得聯(lián)合維數(shù)減少并在基礎(chǔ)RNA表達數(shù)據(jù)上聚類番挺。

DefaultAssay(brain.merge) <- "SCT"
VariableFeatures(brain.merge) <- c(VariableFeatures(brain), VariableFeatures(brain2))
brain.merge <- RunPCA(brain.merge, verbose = FALSE)
brain.merge <- FindNeighbors(brain.merge, dims = 1:30)
brain.merge <- FindClusters(brain.merge, verbose = FALSE)
brain.merge <- RunUMAP(brain.merge, dims = 1:30)

最后，可以在單個UMAP圖中共同可視化數(shù)據(jù)吗浩。SpatialDimPlot()并SpatialFeaturePlot()`默認將所有切片繪制為列，將分組/功能繪制為行没隘。

DimPlot(brain.merge, reduction = "umap", group.by = c("ident", "orig.ident"))

image

SpatialDimPlot(brain.merge)

image

SpatialFeaturePlot(brain.merge, features = c("Hpca", "Plp1"))

image.png

致謝

我們要感謝Nigel Delaney和Stephen Williams對新的空間Seurat代碼的有益反饋和貢獻懂扼。

幻燈片序列

數(shù)據(jù)集

在這里，我們將分析使用小鼠海馬體的Slide-seq v2生成的數(shù)據(jù)集右蒲。本教程將采用與10倍Visium數(shù)據(jù)的空間暈影大致相同的結(jié)構(gòu)阀湿，但專門針對特定于Slide-seq數(shù)據(jù)的演示進行了量身定制。

您可以使用我們的SeuratData包來輕松訪問數(shù)據(jù)瑰妄，如下所示陷嘴。安裝數(shù)據(jù)集后，您可以鍵入?ssHippo以查看用于創(chuàng)建Seurat對象的命令间坐。

InstallData("ssHippo")

slide.seq <- LoadData("ssHippo")

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

通過基因表達數(shù)據(jù)對珠子進行的初始預處理步驟與其他空間Seurat分析和典型的scRNA-seq實驗相似灾挨。在這里，我們注意到許多珠子的UMI計數(shù)特別低竹宋，但選擇保留所有檢測到的珠子用于下游分析劳澄。

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

然后，我們使用sctransform標準化數(shù)據(jù)并執(zhí)行標準的scRNA-seq維數(shù)減少和聚類工作流蜈七。

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空間（使用DimPlot()）或在珠坐標空間使用來可視化聚類的結(jié)果[SpatialDimPlot()](https://satijalab.org/seurat/reference/SpatialPlot.html)。

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

[圖片上傳失敗...(image-30c441-1615977780709)]

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

[圖片上傳失敗...(image-bea41c-1615977780709)]

與scRNA-seq參考整合

為了促進Slide-seq數(shù)據(jù)集的細胞類型注釋飒硅，我們利用了由Saunders *砂缩，Macosko *等人制作的現(xiàn)有小鼠單細胞RNA-seq海馬數(shù)據(jù)集。2018三娩。數(shù)據(jù)可在此處作為已處理的Seurat對象下載庵芭，而原始計數(shù)矩陣可在DropViz網(wǎng)站上獲得。

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

本文的原始注釋在Seurat對象的單元元數(shù)據(jù)中提供雀监。這些注釋以幾種“解決方案”提供喳挑，從大類（ref$class）到細胞類型（）中的子類ref$subcluster。出于本小插圖的目的滔悉，我們將對細胞類型注釋（ref$celltype）進行修改伊诵，使之達到良好的平衡。

我們將從運行Seurat標簽轉(zhuǎn)移方法開始回官，以預測每個珠子的主要細胞類型曹宴。

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

然后，我們可以可視化一些主要預期類別的預測分數(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))

image.png

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

識別空間可變特征

正如Visium插圖中提到的那樣笛坦，我們可以通過兩種通用方法來識別空間可變的特征：預先標注的解剖區(qū)域之間的差異表達測試或測量特征對空間位置的依賴性的統(tǒng)計信息区转。

在這里，我們FindSpatiallyVariableFeatures()通過設(shè)置來實現(xiàn)Moran's I的實現(xiàn)版扩，以演示后者method = 'moransi'废离。Moran's I計算總體空間自相關(guān)，并給出一個統(tǒng)計量（類似于相關(guān)系數(shù)）礁芦，該統(tǒng)計量用于衡量要素對空間位置的依賴性蜻韭。這使我們能夠根據(jù)表達的空間變化程度對要素進行排名。為了便于快速估算此統(tǒng)計信息柿扣，我們實施了一種基本的分級策略肖方，該策略將在Slide-seq圓盤上繪制一個矩形網(wǎng)格，并對每個分級中的特征和位置取平均未状。x和y方向上的倉數(shù)分別由x.cuts和y.cuts參數(shù)控制俯画。此外，雖然不是必需的司草，但安裝可選Rfast2軟件包（install.packages('Rfast2')`）艰垂，將通過更有效的實施方式來顯著減少運行時間。

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)

現(xiàn)在埋虹，我們可視化由Moran's I識別的前6個特征的表達材泄。

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