Seurat4.0系列教程告一段落，但這決不是終點(diǎn)。這個系列教程是給大家打開一扇窗，讓大家知道Seurat4.0有這些功能可用，少走彎路。后續(xù)還要自己深入研究。學(xué)，然后知不足。加油吧！少年！不足之處，懇請大家批評指正！
概述
本教程演示了如何使用 Seurat 4.0 分析單細(xì)胞空間轉(zhuǎn)錄組數(shù)據(jù)。雖然分析流程類似于單細(xì)胞RNA-seq分析的Seurat工作流程，但我們引入了更新的交互和可視化工具，特別強(qiáng)調(diào)了空間和分子信息的整合。本教程將涵蓋以下任務(wù)：
- 標(biāo)準(zhǔn)化
- 降維和聚類
- 檢測空間高變基因
- 交互式可視化
- 與單細(xì)胞RNA-seq數(shù)據(jù)整合
- 處理多個切片
對于第一個教程，我們分析使用Visium 技術(shù)從 10 x genomics生成的數(shù)據(jù)集。我們將在不遠(yuǎn)的將來擴(kuò)展seurat，以處理其他數(shù)據(jù)類型，包括SLIDE-Seq、STARmap和MERFISH。
首先，我們加載Seurat和所需的其他包。
library(Seurat)
library(SeuratData)
library(ggplot2)
library(patchwork)
library(dplyr)
10x Visium
數(shù)據(jù)
這里，我們將使用最近發(fā)布的使用Visium v1 chemistry產(chǎn)生的小鼠大腦切片的數(shù)據(jù)集。
可以在此處下載數(shù)據(jù)，然后使用Load10X_Spatial()將其加載到 Seurat 中。這一步讀取spaceranger流程的輸出，并返回一個包含 spot-level表達(dá)數(shù)據(jù)以及組織切片相關(guān)圖像的 Seurat 對象?？梢允褂?a target="_blank">SeuratData 包輕松訪問數(shù)據(jù)，如下所示。
InstallData("stxBrain")
brain <- LoadData("stxBrain", type = "anterior1")
數(shù)據(jù)預(yù)處理
通過基因表達(dá)數(shù)據(jù)執(zhí)行的預(yù)處理步驟類似于傳統(tǒng)的 scRNA-seq 。首先需要使數(shù)據(jù)標(biāo)準(zhǔn)化，以便較少數(shù)據(jù)之間測序深度的差異。分子計(jì)數(shù)/點(diǎn)的差異對于空間數(shù)據(jù)集來說可能很大，特別是如果整個組織的細(xì)胞密度存在差異。在這里看到巨大的異質(zhì)性，這需要有效的標(biāo)準(zhǔn)化。
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)
這些圖表明，跨點(diǎn)分子計(jì)數(shù)的差異不僅在技術(shù)上，而且依賴于組織解剖學(xué)。例如，神經(jīng)元（如皮質(zhì)白質(zhì)）缺失的組織區(qū)域，可反復(fù)顯示較低的分子計(jì)數(shù)。因此，標(biāo)準(zhǔn)化方法（如函數(shù)LogNormalize()）可能會有問題，這些方法使每個數(shù)據(jù)點(diǎn)在標(biāo)準(zhǔn)化后具有相同的基礎(chǔ)"大小"。
作為替代方案，我們建議使用 sctransform，該模型構(gòu)建了基因表達(dá)的常規(guī)負(fù)二元模型，以便在保持生物差異的同時考慮技術(shù)。sctransform 使數(shù)據(jù)標(biāo)準(zhǔn)化，檢測高變異基因，并將數(shù)據(jù)存儲在SCTassay.
brain <- SCTransform(brain, assay = "Spatial", verbose = FALSE)
基因表達(dá)可視化
在 Seurat 中，我們具有探索空間數(shù)據(jù)固有的視覺性質(zhì)并與之交互的功能。SpatialFeaturePlot()是Seurat 的FeaturePlot()擴(kuò)展，并且可以在組織學(xué)上覆蓋分子數(shù)據(jù)。例如，在小鼠大腦的這一數(shù)據(jù)集中，Hpca基因是一個強(qiáng)大的海馬標(biāo)記，Ttr是choroid plexus的標(biāo)記。
SpatialFeaturePlot(brain, features = c("Hpca", "Ttr"))
Seurat 中的默認(rèn)參數(shù)強(qiáng)調(diào)分子數(shù)據(jù)的可視化。但是，您也可以通過更改以下參數(shù)來調(diào)整點(diǎn)的大小（及其透明度），以改善組織學(xué)圖像的可視化：
-
pt.size.factor-這將擴(kuò)展點(diǎn)的大小。默認(rèn)值為1.6 -
alpha-最低和最大的透明度。默認(rèn)值為 c（1，1）。 - 嘗試設(shè)置為 c（0.1， 1），以降低低表達(dá)點(diǎn)的透明度
alpha
p1 <- SpatialFeaturePlot(brain, features = "Ttr", pt.size.factor = 1)
p2 <- SpatialFeaturePlot(brain, features = "Ttr", alpha = c(0.1, 1))
p1 + p2
降維、聚類和可視化
然后，我們可以繼續(xù)運(yùn)行 RNA 表達(dá)數(shù)據(jù)上的降維和聚類，使用與用于 scRNA-seq 分析相同的工作流。
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)
然后，我們可以使用DimPlot()在 UMAP 空間中或使用SpatialDimPlot()可視化聚類的結(jié)果。
p1 <- DimPlot(brain, reduction = "umap", label = TRUE)
p2 <- SpatialDimPlot(brain, label = TRUE, label.size = 3)
p1 + p2
由于有許多顏色，因此可視化哪個 voxel 屬于哪個聚類可能具有挑戰(zhàn)性。我們有幾個策略來幫助解決此問題。設(shè)置label參數(shù)時，每個集群的中位數(shù)放置一個彩色框（參見上圖）。
您也可以使用cells.highlight參數(shù)在SpatialDimPlot() 。這對于區(qū)分單個群的局部空間非常有用，如下所示：
SpatialDimPlot(brain, cells.highlight = CellsByIdentities(object = brain, idents = c(2, 1, 4, 3,
5, 8)), facet.highlight = TRUE, ncol = 3)
交互式繪圖
我們還建立了一些交互式繪圖功能。SpatialDimPlot()和SpatialDimPlot()兩者現(xiàn)在都有一個參數(shù)，當(dāng)設(shè)置interactive為TRUE，將打開Rstudio窗格與shiny互動。下面的示例演示了一個交互式，您可以在點(diǎn)上懸停并查看細(xì)胞名稱和當(dāng)前marker。
SpatialDimPlot(brain, interactive = TRUE)
對于SpatialFeaturePlot() ，設(shè)置交互式為TRUE，顯示一個交互式窗格，您可以在其中調(diào)整點(diǎn)的透明度、點(diǎn)大小以及正在繪制的函數(shù)。在瀏覽數(shù)據(jù)后，選擇已完成的按鈕將返回最后一個活動圖作為 ggplot 對象。
SpatialFeaturePlot(brain, features = "Ttr", interactive = TRUE)
LinkedDimPlot()函數(shù)將 UMAP 與組織圖像連接，并允許進(jìn)行交互式選擇。例如，您可以選擇 UMAP 繪圖中的區(qū)域，并將突出顯示圖像表示中的相應(yīng)點(diǎn)。
LinkedDimPlot(brain)
空間高變基因的識別
Seurat 提供兩個工作流程來識別與組織內(nèi)空間位置相關(guān)的分子特征。首先是根據(jù)組織內(nèi)預(yù)先注釋的解剖區(qū)域進(jìn)行差異表達(dá)分析，這些區(qū)域可能由不受監(jiān)督的聚類或先驗(yàn)的知識確定。在這種情況下，此策略將起作用，因?yàn)樯鲜黾罕憩F(xià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)
另一種方法，是FindSpatiallyVariables()，是尋找空間模式中沒有預(yù)先注釋的基因。默認(rèn)方法是method = 'markvariogram`。
brain <- FindSpatiallyVariableFeatures(brain, assay = "SCT", features = VariableFeatures(brain)[1:1000],
selection.method = "markvariogram")
現(xiàn)在，我們可視化前 6 個基因的表達(dá)。
top.features <- head(SpatiallyVariableFeatures(brain, selection.method = "markvariogram"), 6)
SpatialFeaturePlot(brain, features = top.features, ncol = 3, alpha = c(0.1, 1))
解剖區(qū)域取子集
與single-cell對象一樣，可以對象取子集。在這里，我們?nèi)〕銮邦~皮層。此過程還有助于將這些數(shù)據(jù)與下一節(jié)中的外皮層 scRNA-seq 數(shù)據(jù)集整合。首先，我們?nèi)∪旱淖蛹?，然后根?jù)確切的位置進(jìn)一步細(xì)分。subset后，我們可以在完整圖像或裁剪圖像上可視化皮質(zhì)細(xì)胞。
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
與single-cell數(shù)據(jù)的整合
在~50um時， visium assay檢測的點(diǎn)將包含多個細(xì)胞的表達(dá)特征。隨著越來越多的可用 scRNA-seq 數(shù)據(jù)，用戶可能有興趣對每個空間 voxels "去卷積"來預(yù)測細(xì)胞類型的潛在組成。在準(zhǔn)備此教程時，我們測試了各種去卷積和整合方法，使用了由 Allen 研究所生成的 14，000 個成年鼠皮質(zhì)細(xì)胞的scRNA-seq 數(shù)據(jù)集。我們發(fā)現(xiàn)使用整合方法（而不是去卷積方法）更優(yōu)越，可能是因?yàn)榭臻g和單細(xì)胞數(shù)據(jù)集的噪聲模型存在較大差異，整合方法針對這些差異進(jìn)行了穩(wěn)健的設(shè)計(jì)。因此，我們應(yīng)用 Seurat v3中引入的基于"錨點(diǎn)"的sctransform 整合工作流，從而能夠?qū)⒓?xì)胞從參考集轉(zhuǎn)到查詢集進(jìn)行注釋。
首先加載數(shù)據(jù)（此處下載），預(yù)處理 scRNA-seq 參考集，然后執(zhí)行標(biāo)簽轉(zhuǎn)移。該過程為每個點(diǎn)輸出每個scRNA-seq細(xì)胞的概率分類。我們將這些預(yù)測添加到新的Seurat 對象進(jìn)行分析。
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)
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)在，我們得到每個類每個點(diǎn)的預(yù)測分?jǐn)?shù)。最感興趣的是frontal cortex region的laminar excitatory neurons。在這里，我們可以區(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)
基于這些預(yù)測分?jǐn)?shù)，我們還可以預(yù)測其位置受空間限制的細(xì)胞類型。使用基于標(biāo)記點(diǎn)過程的相同方法來定義空間變異基因，但使用細(xì)胞類型預(yù)測分?jǐn)?shù)作為"標(biāo)記"，而不是基因表達(dá)。
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)
最后，我們表明，我們的整合程序能夠恢復(fù)神經(jīng)元和非神經(jīng)元子集已知的空間定位模式，包括laminar excitatory, layer-1 astrocytes, and the cortical grey matter。
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))
在seurat處理多個切片
鼠大腦的這個數(shù)據(jù)集包含另一個與大腦另一半相對應(yīng)的切片。在這里，我們讀入它，并執(zhí)行相同的初始標(biāo)準(zhǔn)化。
brain2 <- LoadData("stxBrain", type = "posterior1")
brain2 <- SCTransform(brain2, assay = "Spatial", verbose = FALSE)
為了處理同一 Seurat 對象中的多個切片，我們提供了merge。
brain.merge <- merge(brain, brain2)
這便能在下面的RNA表達(dá)數(shù)據(jù)上實(shí)現(xiàn)整合降維和聚類。
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)
最后，數(shù)據(jù)可以在單個 UMAP 繪圖中共同可視化。 默認(rèn)情況下，將所有切片繪制為列，分組/基因作為行。
DimPlot(brain.merge, reduction = "umap", group.by = c("ident", "orig.ident"))
SpatialDimPlot(brain.merge)
SpatialFeaturePlot(brain.merge, features = c("Hpca", "Plp1"))
Slide-seq
數(shù)據(jù)
這里，我們將分析使用小鼠海馬的Slide-seq v2數(shù)據(jù)集。此教程將遵循與 10x Visium 數(shù)據(jù)的空間教程大致相同的結(jié)構(gòu)。
您可以使用SeuratData 包輕松訪問數(shù)據(jù)，如下所示。安裝數(shù)據(jù)集后，您可以鍵入?ssHippo以查看用于創(chuàng)建 Seurat 對象的命令。
InstallData("ssHippo")
slide.seq <- LoadData("ssHippo")
數(shù)據(jù)預(yù)處理
通過基因表達(dá)數(shù)據(jù)對bead 的初始預(yù)處理步驟類似于其他空間seurat分析和傳統(tǒng)的scRNA-seq實(shí)驗(yàn)。在這里，我們注意到許多beads 包含特別低的 UMI 計(jì)數(shù)，我們選擇保留所有檢測到的beads 進(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)
然后，我們使用sctransform進(jìn)行數(shù)據(jù)標(biāo)準(zhǔn)化，并執(zhí)行標(biāo)準(zhǔn)的 scRNA-seq 降維和聚類工作流。
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 space或bead coordinate space中可視化聚類的結(jié)果。
plot1 <- DimPlot(slide.seq, reduction = "umap", label = TRUE)
plot2 <- SpatialDimPlot(slide.seq, stroke = 0)
plot1 + plot2
SpatialDimPlot(slide.seq, cells.highlight = CellsByIdentities(object = slide.seq, idents = c(1,
6, 13)), facet.highlight = TRUE)
與scRNA-seq參考整合
為了促進(jìn)Slide-seq數(shù)據(jù)集的細(xì)胞類型注釋，我們利用了由Saunders *，Macosko *等人制作的現(xiàn)有小鼠單細(xì)胞RNA-seq海馬數(shù)據(jù)集。數(shù)據(jù)可在此處作為已處理的Seurat對象下載，而原始計(jì)數(shù)矩陣可在DropViz網(wǎng)站上獲得。
ref <- readRDS("../data/mouse_hippocampus_reference.rds")
本文的原始注釋在Seurat對象的metadata 中提供。這些注釋以幾種“resolutions”提供，從大類（ref$class）到亞群ref$subcluster。出于本教程的目的，我們將對細(xì)胞類型注釋（ref$celltype）進(jìn)行修改，使之達(dá)到良好的平衡。
我們將從運(yùn)行Seurat標(biāo)簽轉(zhuǎn)移方法開始，以預(yù)測每個bead的主要細(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)
識別空間變異基因
正如Visium教程中提到的那樣，我們可以通過兩種通用方法來識別空間變異基因：預(yù)先標(biāo)注的解剖區(qū)域之間的差異表達(dá)檢測或基因特征對空間位置的依賴性的統(tǒng)計(jì)信息。
在這里，我們演示后者，通過在FindSpatiallyVariableFeatures()中設(shè)置method = 'moransi'。 Moran's I計(jì)算總體空間自相關(guān)性，并給出一個統(tǒng)計(jì)量（類似于相關(guān)系數(shù)），該統(tǒng)計(jì)量用于衡量要素對空間位置的獨(dú)立性。這使我們能夠根據(jù)基因表達(dá)的變化對基因進(jìn)行排名。為了便于快速估算此統(tǒng)計(jì)信息，我們實(shí)施了一種基本的分級策略，該策略將在Slide-seqpuck 上繪制一個矩形網(wǎng)格，并對每個分級中的基因和位置取平均值。x和y方向上的bins 分別由x.cuts和y.cuts參數(shù)控制。此外，雖然不是必需的，但安裝可選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)
現(xiàn)在，我們可視化由Moran's I識別的前6個基因的表達(dá)。
SpatialFeaturePlot(slide.seq, features = head(SpatiallyVariableFeatures(slide.seq, selection.method = "moransi"),
6), ncol = 3, alpha = c(0.1, 1), max.cutoff = "q95")
