0. 原理簡介
空間轉(zhuǎn)錄組的概念：將基因表達與組織切片的免疫組化圖像相結(jié)合尔破，從而將組織內(nèi)不同細胞的基因表達信息定位到組織的原始空間上，達到直觀檢測組織中不同部位基因表達差異的目的〖崴祝可以避免單細胞轉(zhuǎn)錄組樣本消化對細胞的影響祸泪，結(jié)合組織定位也可以更好的判斷細胞類型。
核心：單細胞轉(zhuǎn)錄組測序+組織定位
空轉(zhuǎn)芯片：
? 一個載玻片可容納4個切片(捕獲區(qū)域)腹忽，制備4個文庫来累；
? 每個6.5mm x 6.5mm捕獲區(qū)域, 有5000個barcode標記的ST位點；
? 每個ST位點有百萬個UMI探針窘奏；
? 每個位點直徑55μm嘹锁，點與點之間相距100μm;
? 每個位點覆蓋1-10個細胞，檢測數(shù)千個基因;
10x公司的Visium空間轉(zhuǎn)錄組方法與10x的單細胞測序在原理上很相似着裹。單細胞測序是用膠珠和油包水方法领猾，把細胞分隔開，同時又用DNA條形碼保留單細胞信息骇扇。Visium空間轉(zhuǎn)錄組則是把切片在芯片上展開摔竿，在空間上用條形碼來保留切片上每個小點的空間位置信息。
1. 載入數(shù)據(jù)
這里使用的是Visium v1得到的小鼠腦組織切片演示數(shù)據(jù)集少孝。包含兩個連續(xù)的anterior切片和2個連續(xù)的posterior切片继低。如果是自己的數(shù)據(jù)，可以使用 Load10X_Spatial() 函數(shù)讀入（filtered_feature_bc_matrix.h5）稍走。
library(Seurat)
library(SeuratData)
library(ggplot2)
library(patchwork)
library(dplyr)
InstallData("stxBrain")
brain <- LoadData("stxBrain", type = "anterior1")
#只讀取了anterior1這個sample的數(shù)據(jù)
我們先來看一下空間轉(zhuǎn)錄組數(shù)據(jù)的Seurat對象
有2696個spots（一共是5000個）
View(brain)
2. 數(shù)據(jù)預(yù)處理
沒有做質(zhì)控哦（過濾方法跟單細胞轉(zhuǎn)錄組一樣的袁翁，過濾掉的細胞在切片上會看到spots變灰了）
brain[['percent.mt']] <- PercentageFeatureSet(brain,pattern = '^mt-')
brain[['percent.rb']] <- PercentageFeatureSet(brain,pattern = '^Rp[ls]')
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測得reads數(shù)的變化不僅是技術(shù)上的問題，而且還取決于組織的解剖結(jié)構(gòu)。As a result, standard approaches (such as the LogNormalize function), which force each data point to have the same underlying ‘size’ after normalization, can be problematic.
作為替代方案，我們建議使用sctransform戚哎，它構(gòu)建了基因表達的正則化負二項式模型，以便在保留生物學(xué)差異的同時考慮技術(shù)偽像掀抹。sctransform可以對數(shù)據(jù)進行歸一化，檢測高變異特征并將數(shù)據(jù)存儲在SCTassay中心俗。
brain <- SCTransform(brain, assay = "Spatial", verbose = FALSE)
拓展：
How do results compare to log-normalization
為了探究標準化方式的差異傲武，作者檢測了sctransform 和log normalization和UMIs的相關(guān)性
# 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
可以看到SCTransform表現(xiàn)的更好
3. 基因表達可視化
SpatialFeaturePlot(brain, features = c("Hpca", "Ttr"))
p1 <- SpatialFeaturePlot(brain, features = "Ttr", pt.size.factor = 1)
p2 <- SpatialFeaturePlot(brain, features = "Ttr", alpha = c(0.1, 1))
p1 + p2
4. 降維蓉驹，聚類和可視化
和scRNA-seq一模一樣（Louvain聚類）
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)
這種聚類方法最容易實現(xiàn)，但是不夠精確揪利。（BayesSpace聚類等效果會更好）
p1 <- DimPlot(brain, reduction = "umap", label = TRUE)
p2 <- SpatialDimPlot(brain, label = TRUE, label.size = 3)
p3 <- SpatialDimPlot(brain, label = TRUE, label.size = 3)&ggsci::scale_fill_igv()
p1|p2|p3 #對接ggsci自定義顏色
SpatialDimPlot(brain, cells.highlight = CellsByIdentities(object = brain, idents = c(2, 1, 4, 3,
5, 8)), facet.highlight = TRUE, ncol = 3)
# 相當于單細胞的split.by
5. 鑒定空間高變基因
FindMarkers找的是上面不同分群間的差異基因
de_markers <- FindMarkers(brain, ident.1 = 5, ident.2 = 6)
View(de_markers)
SpatialFeaturePlot(object = brain, features = rownames(de_markers)[1:3], alpha = c(0.1, 1), ncol = 3)
FindSpatiallyVariableFeatures找的是空間高變基因
空間高變基因指的是表達與空間位置相關(guān)的基因态兴，也可以說是具有空間偏好性的基因∨蔽唬空間轉(zhuǎn)錄組學(xué)允許研究人員調(diào)查基因表達趨勢如何在空間上變化瞻润，從而確定基因表達的空間模式。尋找空間高變基因的算法有很多甜刻，比如SPARK绍撞，SpatialDE等。
brain <- FindSpatiallyVariableFeatures(brain, assay = "SCT", features = VariableFeatures(brain)[1:1000],
selection.method = "markvariogram")
top.features <- head(SpatiallyVariableFeatures(brain, selection.method = "markvariogram"), 6)
SpatialFeaturePlot(brain, features = top.features, ncol = 3, alpha = c(0.1, 1))
6. 細分解剖區(qū)域
這個其實用Loupe Browser手動圈會更方便一點
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
7. 和單細胞數(shù)據(jù)整合
導(dǎo)入單細胞數(shù)據(jù)
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)
對空間數(shù)據(jù)里面選出來的亞群（brain里面挑出來的cortex）數(shù)據(jù)重新進行SCT標準化
# 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)
整合單細胞和空間組學(xué)數(shù)據(jù)（基于錨點轉(zhuǎn)換的去卷積）
注意這里是按組織結(jié)構(gòu)去做注釋的（從空轉(zhuǎn)數(shù)據(jù)中提出了cortex）得院。一般不建議使用大的單細胞數(shù)據(jù)集去做參考傻铣。原因一是速度會很慢，二是容易引入批次效應(yīng)讓結(jié)果不準確祥绞。所以盡量還是找一對一匹配的數(shù)據(jù)集做錨點轉(zhuǎn)換的去卷積非洲，也就是測了空轉(zhuǎn)的數(shù)據(jù)，最好有相應(yīng)的同一個樣本的單細胞數(shù)據(jù)做映射就谜。整合了多個樣品的單細胞數(shù)據(jù)集雖然可能細胞類型更全，但會引入批次效應(yīng)里覆。
##尋找錨點
anchors <- FindTransferAnchors(reference = allen_reference, query = cortex, normalization.method = "SCT")
##錨點轉(zhuǎn)換
predictions.assay <- TransferData(anchorset = anchors, refdata = allen_reference$subclass, prediction.assay = TRUE,
weight.reduction = cortex[["pca"]], dims = 1:30)
cortex[["predictions"]] <- predictions.assay
DefaultAssay(cortex) <- "predictions"
SpatialFeaturePlot(cortex, features = c("L2/3 IT", "L4"), pt.size.factor = 1.6, ncol = 2, crop = TRUE)
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)
SpatialFeaturePlot(cortex, features = c("Astro", "L2/3 IT", "L4", "L5 PT", "L5 IT", "L6 CT", "L6 IT",
"L6b", "Oligo"), pt.size.factor = 1, ncol = 5, crop = FALSE, alpha = c(0.1, 1))
可以看到不同的單細胞細胞群在空間數(shù)據(jù)上的映射都存在一個評分喧枷。這個評分其實不是可信度虹统，可以將其理解為細胞占比。我們知道每個spot其實是多個細胞的混合隧甚，一個spot中有些細胞占比多车荔，有些細胞占比少。比如說上圖第一行第三張圖的L4范圍只有0.1-0.3戚扳，就說明最紅的這些spot中忧便，每個spot中可能最多也只有30%的細胞是L4，其他70%可能是其他細胞帽借，需要結(jié)合其他細胞群的圖來一起看珠增。
參考：https://satijalab.org/seurat/articles/spatial_vignette.html
