思考：來自運(yùn)來大佬的筆記(現(xiàn)在走的路都是大佬為我們鋪好的路，要學(xué)會(huì)珍惜鸭限！)
- 如何將空間數(shù)據(jù)與表達(dá)數(shù)據(jù)關(guān)聯(lián)在一起蜕径？
- 有了空間轉(zhuǎn)錄組數(shù)據(jù)，如何與單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)聯(lián)用里覆？
- 做了多層切片如何展示真實(shí)的三維空間的轉(zhuǎn)錄本信息丧荐？
本教程演示了如何使用Seurat(>=3.2)來分析空間解析的RNA-seq數(shù)據(jù)。雖然分析流程類似于用于單細(xì)胞RNA-seq分析的Seurat工作流喧枷，但我們引入了更新的交互和可視化工具，特別強(qiáng)調(diào)空間和分子信息的集成弓坞。本教程將涵蓋以下任務(wù)隧甚，我們相信這將是常見的許多空間分析:
- 歸一化
- 降維和聚類
- 檢測spatially-variable特性
- 交互式可視化
- 與單細(xì)胞RNA-seq數(shù)據(jù)集成
- 處理多個(gè)切片
對(duì)于我們的第一個(gè)教程，我們分析了與Visium技術(shù)從10x基因組產(chǎn)生的一個(gè)數(shù)據(jù)集渡冻。在不久的將來戚扳，我們將擴(kuò)展Seurat以處理更多的數(shù)據(jù)類型，包括 SLIDE-Seq, STARmap, 和MERFISH.
1.加載包
library(Seurat)
library(SeuratData)
library(ggplot2)
library(patchwork)
library(dplyr)
2.數(shù)據(jù)集
在這里族吻，我們將使用最近發(fā)布的使用Visium v1化學(xué)試劑生成的數(shù)據(jù)集矢狀面小鼠腦切片帽借。有兩個(gè)連續(xù)的前部和兩個(gè)(匹配的)連續(xù)的后部。
您可以從這里下載數(shù)據(jù)超歌，并使用Load10X_Spatial函數(shù)將其加載到Seurat砍艾。這將讀入spaceranger流程的輸出，并返回一個(gè)Seurat對(duì)象巍举，該對(duì)象包含位置級(jí)表達(dá)式數(shù)據(jù)以及組織切片的關(guān)聯(lián)圖像脆荷。您還可以使用我們的SeuratData包來方便地進(jìn)行數(shù)據(jù)訪問，如下所示。安裝數(shù)據(jù)集之后蜓谋，您可以輸入?stxBrain以了解更多信息
InstallData("stxBrain")
brain <- LoadData("stxBrain", type = "anterior1")
brain
An object of class Seurat
31053 features across 2696 samples within 1 assay
Active assay: Spatial (31053 features, 0 variable features)
?
brain@assays$Spatial[1:4,1:4]
AAACAAGTATCTCCCA-1 AAACACCAATAACTGC-1 AAACAGAGCGACTCCT-1 AAACAGCTTTCAGAAG-1
Xkr4 . . . .
Gm1992 . . . .
Gm37381 . . . .
Rp1 . . . .
數(shù)據(jù)包含了31053個(gè)基因梦皮，2696個(gè)spot
來自10x的visium數(shù)據(jù)包括以下數(shù)據(jù)類型，數(shù)據(jù)存儲(chǔ)結(jié)構(gòu)如下：
- 一個(gè)點(diǎn)基因表達(dá)矩陣：每一行為一個(gè)基因桃焕，每一列為一個(gè)spot剑肯，每個(gè)spot可能包含多個(gè)細(xì)胞。
- 組織切片圖像(數(shù)據(jù)采集時(shí)通過H&E染色獲得)
- 將原始高分辨率圖像與此處用于可視化的低分辨率圖像關(guān)聯(lián)的縮放因子观堂。
在這里spot表達(dá)矩陣是需要重點(diǎn)理解的：
在Seurat實(shí)驗(yàn)對(duì)象中退子，spot基因表達(dá)矩陣類似于典型的“RNA”實(shí)驗(yàn)，但包含了spot水平型将，而不是單細(xì)胞水平的數(shù)據(jù)寂祥。映像本身存儲(chǔ)在Seurat對(duì)象中的一個(gè)新的映像槽中。圖像插槽還存儲(chǔ)了必要的信息七兜，以關(guān)聯(lián)的位置與其在組織圖像上的物理位置丸凭。
3.數(shù)據(jù)預(yù)處理
spot水平的基因表達(dá)數(shù)據(jù)的初始預(yù)處理步驟與典型的scRNA-seq實(shí)驗(yàn)相似。我們首先需要對(duì)數(shù)據(jù)進(jìn)行歸一化腕铸，以解釋數(shù)據(jù)點(diǎn)間序列深度的差異惜犀。我們注意到，分子counts/spot的差異對(duì)于空間數(shù)據(jù)集來說是實(shí)質(zhì)性的狠裹，特別是當(dāng)組織中的細(xì)胞密度存在差異時(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)
這些圖表明莉御，不同spot的分子計(jì)數(shù)的差異不僅是技術(shù)上的，而且還取決于組織解剖學(xué)俗冻。例如礁叔，神經(jīng)元減少的組織區(qū)域(如皮質(zhì)白質(zhì)cortical white matter)，重復(fù)地表現(xiàn)出較低的分子計(jì)數(shù)迄薄。因此琅关，標(biāo)準(zhǔn)方法(如LogNormalize函數(shù))在標(biāo)準(zhǔn)化后強(qiáng)制每個(gè)數(shù)據(jù)點(diǎn)具有相同的底層“大小”，這可能會(huì)帶來問題讥蔽。
作為一種替代方法涣易，我們建議使用sctransform (Hafemeister和Satija, Genome Biology 2019)，該方法構(gòu)建正則化的負(fù)二項(xiàng)基因表達(dá)模型冶伞，以解釋技術(shù)人工產(chǎn)物新症，同時(shí)保留生物方差。有關(guān)sctransform的更多細(xì)節(jié)碰缔，請(qǐng)參見這里的論文和這里的Seurat教程账劲。sctransform對(duì)數(shù)據(jù)進(jìn)行規(guī)范化，檢測高方差特征，并將數(shù)據(jù)存儲(chǔ)在SCT assay中瀑焦。
brain <- SCTransform(brain, assay = "Spatial", verbose = FALSE)
注意這里有一個(gè)默認(rèn)參數(shù)：variable.features.n = 3000腌且，選擇的是3000個(gè)基因。
SCTransform與LogNormalize兩種標(biāo)化方法相比：
為了探究歸一化方法的差異榛瓮，我們研究了sctransform和log歸一化結(jié)果如何與UMIs的數(shù)量相關(guān)聯(lián)铺董。為了進(jìn)行比較，我們首先重新運(yùn)行sctransform來存儲(chǔ)所有基因的值禀晓，并通過NormalizeData運(yùn)行LogNormalize過程精续。
# rerun normalization to store sctransform residuals for all genes
brain <- SCTransform(brain, assay = "Spatial", return.only.var.genes = FALSE, verbose = FALSE)
# also run standard log normalization for comparison
brain <- NormalizeData(brain, verbose = FALSE, assay = "Spatial")
?
# Computes the correlation of the log normalized data and sctransform residuals with the number
# of UMIs
brain <- GroupCorrelation(brain, group.assay = "Spatial", assay = "Spatial", slot = "data", do.plot = FALSE)
brain <- GroupCorrelation(brain, group.assay = "Spatial", assay = "SCT", slot = "scale.data", do.plot = FALSE)
?
p1 <- GroupCorrelationPlot(brain, assay = "Spatial", cor = "nCount_Spatial_cor") + ggtitle("Log Normalization") +
theme(plot.title = element_text(hjust = 0.5))
p2 <- GroupCorrelationPlot(brain, assay = "SCT", cor = "nCount_Spatial_cor") + ggtitle("SCTransform Normalization") +
theme(plot.title = element_text(hjust = 0.5))
p1 + p2
對(duì)于上面的箱線圖，我們計(jì)算每個(gè)特征(基因)與UMIs(這里的nCount_Spatial變量)數(shù)量的相關(guān)性粹懒。然后重付，我們根據(jù)基因的平均表達(dá)水平將其分組，并生成這些相關(guān)性的箱線圖凫乖。你可以看到确垫，在前三組LogNormalize未能充分正常化基因帽芽，這表明技術(shù)因素繼續(xù)影響標(biāo)準(zhǔn)化表達(dá)估計(jì)高表達(dá)基因删掀。相反，sctransform標(biāo)準(zhǔn)化實(shí)質(zhì)上減輕了這種影響导街。
4.基因表達(dá)的可視化
在Seurat v3.2中披泪，我們添加了一些新功能，用于探索空間數(shù)據(jù)固有的視覺特性并與之交互搬瑰。Seurat的空間特征圖功能擴(kuò)展了特征圖款票，可以將分子數(shù)據(jù)覆蓋在組織組織學(xué)上。例如跌捆，在這組小鼠大腦的數(shù)據(jù)中徽职，Hpca基因是一個(gè)強(qiáng)大的海馬標(biāo)記，而Ttr是脈絡(luò)膜叢的標(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：因子-這將縮放斑點(diǎn)的大小。默認(rèn)是1.6
- alpha：最小和最大透明度陶冷。默認(rèn)為c(1,1)
嘗試設(shè)置為alpha c(0.1, 1)钙姊，來降低表達(dá)式較低的點(diǎn)的透明度
p1 <- SpatialFeaturePlot(brain, features = "Ttr", pt.size.factor = 1)
p2 <- SpatialFeaturePlot(brain, features = "Ttr", alpha = c(0.1, 1))
p1 + p2
5.降維、聚類和可視化
然后埂伦，我們可以繼續(xù)對(duì)RNA表達(dá)數(shù)據(jù)進(jìn)行降維和聚類煞额，使用與我們用于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)
然后，我們可以在UMAP空間(使用DimPlot)或在圖像上使用SpatialDimPlot疊加顯示聚類結(jié)果
p1 <- DimPlot(brain, reduction = "umap", label = TRUE)
p2 <- SpatialDimPlot(brain, label = TRUE, label.size = 3)
p1 + p2
由于有許多顏色膊毁，它可能是具有挑戰(zhàn)性的可視化哪個(gè)體素屬于哪個(gè)cluster胀莹。我們有一些策略來幫助解決這個(gè)問題。設(shè)置label參數(shù)會(huì)在每個(gè)集群的中間位置放置一個(gè)彩色框(見上面的圖)婚温。
您還可以使用cells.highlight參數(shù)來標(biāo)定空間圖上感興趣的特定單元格描焰。這對(duì)于區(qū)分單個(gè)集群的空間定位非常有用，如下所示：
SpatialDimPlot(brain, cells.highlight = CellsByIdentities(object = brain, idents = c(2, 1, 4, 3, 5, 8)), facet.highlight = TRUE, ncol = 3)
6.識(shí)別空間變量特征
Seurat提供了兩個(gè)工作流程來識(shí)別組織內(nèi)與空間位置相關(guān)的分子特征栅螟。第一種是基于組織內(nèi)預(yù)先注釋的解剖區(qū)域來進(jìn)行差異表達(dá)荆秦，這可以通過無監(jiān)督聚類或先驗(yàn)知識(shí)來確定。在這種情況下力图，這個(gè)策略是有效的步绸，因?yàn)樯厦娴腸luster有明確的空間限制。
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中實(shí)現(xiàn)的另一種方法是在沒有預(yù)先注釋的情況下搜索顯示空間模式的特性吃媒。默認(rèn)的方法(method = 'markvariogram ')受到 Trendsceek的啟發(fā)瓤介，后者將空間轉(zhuǎn)錄組數(shù)據(jù)建模為一個(gè)標(biāo)記點(diǎn)過程，并計(jì)算一個(gè)'variogram'晓折，它識(shí)別基因的表達(dá)水平取決于它們的空間位置惑朦。更具體地說，這個(gè)過程計(jì)算gamma(r)值漓概，測量的是距離為一定“r”的兩點(diǎn)之間的依賴關(guān)系漾月。默認(rèn)情況下，我們在這些分析中使用的r值為“5”胃珍，并且僅為可變基因計(jì)算這些值(其中變異的計(jì)算獨(dú)立于空間位置)以節(jié)省時(shí)間梁肿。
我們注意到在文獻(xiàn)中有多種方法來完成這項(xiàng)任務(wù)，包括空間法SpatialDE和斑點(diǎn)法Splotch觅彰。我們鼓勵(lì)有興趣的用戶探索這些方法吩蔑，并希望在不久的將來增加對(duì)它們的支持。
brain <- FindSpatiallyVariableFeatures(brain, assay = "SCT", features = VariableFeatures(brain)[1:1000], selection.method = "markvariogram")
現(xiàn)在我們將通過這種方法確定的前6個(gè)特征的表達(dá)形象化
top.features <- head(SpatiallyVariableFeatures(brain, selection.method ="markvariogram"), 6)
SpatialFeaturePlot(brain, features = top.features, ncol = 3, alpha = c(0.1, 1))
7.劃分解剖區(qū)域
與單細(xì)胞對(duì)象一樣填抬，可以對(duì)對(duì)象進(jìn)行子集化烛芬，以集中于數(shù)據(jù)的子集。這里飒责，我們大致劃分了額葉皮層赘娄。這個(gè)過程還將在下一節(jié)中促進(jìn)這些數(shù)據(jù)與皮質(zhì)scRNA-seq數(shù)據(jù)集的集成。首先選取cluster的子集宏蛉，然后根據(jù)具體位置進(jìn)一步細(xì)分遣臼。經(jīng)過亞設(shè)置硼一，我們可以在完整的圖像或裁剪的圖像上看到皮質(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
8.與單細(xì)胞數(shù)據(jù)整合
在50um時(shí)醉旦，來自visium的斑點(diǎn)將包圍多個(gè)細(xì)胞的表達(dá)譜。對(duì)于不斷增長的可用scRNA-seq數(shù)據(jù)的系統(tǒng)列表肥矢，用戶可能會(huì)對(duì)“反褶積”每個(gè)空間體素感興趣族阅，以預(yù)測細(xì)胞類型的底層組成认臊。在準(zhǔn)備這個(gè)教程時(shí)扳炬，我們測試了各種各樣的decovonlution和集成方法，使用參考的reference scRNA-seq dataset西采，包括來自Allen研究所的約14,000成年小鼠皮質(zhì)細(xì)胞分類凰萨，由SMART-Seq2生成。我們一直發(fā)現(xiàn)使用積分方法(相對(duì)于反褶積方法)會(huì)有更好的性能械馆，這可能是因?yàn)槊枋隹臻g和單細(xì)胞數(shù)據(jù)集的噪聲模型存在本質(zhì)上的差異胖眷，而積分方法是專門設(shè)計(jì)來應(yīng)對(duì)這些差異的。因此我們應(yīng)用“錨”的集成工作流最近在修v3,介紹,使注釋的概率轉(zhuǎn)移從一個(gè)參考查詢集霹崎。因此,我們遵循這里介紹的標(biāo)簽轉(zhuǎn)移流程,利用sctransform正成翰螅化,但預(yù)測新方法被開發(fā)來完成這項(xiàng)任務(wù)。
我們首先加載數(shù)據(jù)(這里可以下載)尾菇，預(yù)處理scRNA-seq參考數(shù)據(jù)集境析，然后執(zhí)行標(biāo)簽傳輸。該過程為每個(gè)點(diǎn)輸出每個(gè)scRNA-seq派生類的概率分類派诬。我們在Seurat對(duì)象中加入這些預(yù)測作為一個(gè)新的分析劳淆。
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"]])
cortex[["predictions"]] <- predictions.assay
現(xiàn)在我們得到每個(gè)class的每個(gè)spot的預(yù)測分?jǐn)?shù)。在額葉皮質(zhì)區(qū)域特別有趣的是層狀興奮神經(jīng)元(laminar excitatory neurons)默赂。在這里沛鸵，我們可以區(qū)分這些神經(jīng)元亞型的不同順序?qū)樱?
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)亞群的已知空間定位模式，包括層興奮性奖恰、第一層星形膠質(zhì)細(xì)胞和皮質(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))
此處圖比較大，去看官網(wǎng)吧瑟啃。
9.在Seurat中使用多個(gè)切片
這個(gè)老鼠大腦的數(shù)據(jù)集包含了另一個(gè)對(duì)應(yīng)于大腦的另一半的切片趾徽。在這里，我們讀取它并執(zhí)行相同的初始化翰守。
brain2 <- LoadData("stxBrain", type = "posterior1")
brain2 <- SCTransform(brain2, assay = "Spatial", verbose = FALSE)
為了在同一個(gè)Seurat對(duì)象中使用多個(gè)片，我們提供了merge函數(shù)疲酌。
brain.merge <- merge(brain, brain2)
這樣就可以減少潛在的RNA表達(dá)數(shù)據(jù)的聯(lián)合維數(shù)和聚類蜡峰。
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ù)可以聯(lián)合顯示在一個(gè)UMAP圖中。SpatialDimPlot和SpatialFeaturePlot默認(rèn)情況下將所有的切片繪制為列湿颅，將分組/特性繪制為行载绿。
DimPlot(brain.merge, reduction = "umap", group.by = c("ident", "orig.ident"))
SpatialDimPlot(brain.merge)
SpatialFeaturePlot(brain.merge, features = c("Hpca", "Plp1"))
10.致謝
我們要感謝Nigel Delaney和Stephen Williams，感謝他們對(duì)新的空間Seurat代碼的反饋和貢獻(xiàn)油航。
