引言
在這篇指南中，我們介紹了Seurat的一個新擴(kuò)展功能，用以分析新型的空間解析數(shù)據(jù)，將重點介紹由不同成像技術(shù)生成的三個公開數(shù)據(jù)集。
- Vizgen MERSCOPE（用于小鼠大腦研究）
- Nanostring CosMx空間分子成像儀（用于FFPE人類肺組織）
- Akoya CODEX（用于人類淋巴結(jié)研究）
小鼠大腦：10x Genomics Xenium In Situ
在本節(jié)中，我們將對Xenium平臺生成的數(shù)據(jù)進(jìn)行分析。本指南展示了如何導(dǎo)入Xenium平臺輸出的每轉(zhuǎn)錄本位置信息、細(xì)胞與基因的表達(dá)矩陣、細(xì)胞邊界分割以及細(xì)胞中心點數(shù)據(jù)。得到的Seurat對象將包含每個細(xì)胞的基因表達(dá)概況、細(xì)胞的中心點和邊界位置，以及每個被檢測到的轉(zhuǎn)錄本的具體位置。這些細(xì)胞層面的基因表達(dá)概況與標(biāo)準(zhǔn)的單細(xì)胞RNA測序數(shù)據(jù)相似，可以運用相同的分析工具進(jìn)行處理。
這里使用的是10x Genomics公司為Xenium Explorer Demo提供的“Tiny”子集數(shù)據(jù)集，該數(shù)據(jù)集源自“Fresh Frozen Mouse Brain”，下載方式如下。這些分析步驟同樣適用于更大的“完整冠狀切面”數(shù)據(jù)集，但處理時間會更長。
wget https://cf.10xgenomics.com/samples/xenium/1.0.2/Xenium_V1_FF_Mouse_Brain_Coronal_Subset_CTX_HP/Xenium_V1_FF_Mouse_Brain_Coronal_Subset_CTX_HP_outs.zip
unzip Xenium_V1_FF_Mouse_Brain_Coronal_Subset_CTX_HP_outs.zip
首先，我們讀取數(shù)據(jù)集并創(chuàng)建一個 Seurat 對象。提供 Xenium 運行的數(shù)據(jù)文件夾的路徑作為輸入路徑。 RNA 數(shù)據(jù)存儲在 Seurat 對象的 Xenium 分析中。此步驟大約需要一分鐘。
path <- "/brahms/hartmana/vignette_data/xenium_tiny_subset"
# Load the Xenium data
xenium.obj <- LoadXenium(path, fov = "fov")
# remove cells with 0 counts
xenium.obj <- subset(xenium.obj, subset = nCount_Xenium > 0)
空間數(shù)據(jù)被嵌入到Seurat對象的特定字段中，這些字段以加載的“視野”（Field of View，簡稱FOV）名稱來標(biāo)識。一開始，所有的數(shù)據(jù)都被加載到默認(rèn)的FOV中，名為fov。接下來，我們將創(chuàng)建一個新的視野，這個視野將聚焦于我們特別關(guān)注的區(qū)域。
Xenium分析工具提供了Seurat的標(biāo)準(zhǔn)質(zhì)量控制圖表。例如，我們可以通過小提琴圖來查看每個細(xì)胞中的基因數(shù)量（nFeature_Xenium）和每個細(xì)胞中的轉(zhuǎn)錄本計數(shù)（nCount_Xenium）。這些圖表幫助我們對數(shù)據(jù)集的質(zhì)量有一個直觀的了解。
VlnPlot(xenium.obj, features = c("nFeature_Xenium", "nCount_Xenium"), ncol = 2, pt.size = 0)
接下來，我們使用 ImageDimPlot() 在組織上繪制全抑制神經(jīng)元標(biāo)記 Gad1、抑制神經(jīng)元亞型標(biāo)記 Pvalb 和 Sst 以及星形膠質(zhì)細(xì)胞標(biāo)記 Gfap 的位置。
ImageDimPlot(xenium.obj, fov = "fov", molecules = c("Gad1", "Sst", "Pvalb", "Gfap"), nmols = 20000)
在這里，我們使用 ImageFeaturePlot() 在每個細(xì)胞水平上可視化一些關(guān)鍵層標(biāo)記基因的表達(dá)水平，這類似于用于可視化 2D 嵌入表達(dá)的 FeaturePlot() 函數(shù)。我們手動將每個基因的 max.cutoff 調(diào)整到其計數(shù)分布的大約第 90 個百分位數(shù)（可以使用 max.cutoff='q90' 指定）以提高對比度。
ImageFeaturePlot(xenium.obj, features = c("Cux2", "Rorb", "Bcl11b", "Foxp2"), max.cutoff = c(25,
35, 12, 10), size = 0.75, cols = c("white", "red"))
我們可以使用 Crop() 函數(shù)放大所選區(qū)域。放大后，我們可以可視化細(xì)胞分割邊界以及單個分子。
cropped.coords <- Crop(xenium.obj[["fov"]], x = c(1200, 2900), y = c(3750, 4550), coords = "plot")
xenium.obj[["zoom"]] <- cropped.coords
# visualize cropped area with cell segmentations & selected molecules
DefaultBoundary(xenium.obj[["zoom"]]) <- "segmentation"
ImageDimPlot(xenium.obj, fov = "zoom", axes = TRUE, border.color = "white", border.size = 0.1, cols = "polychrome",
coord.fixed = FALSE, molecules = c("Gad1", "Sst", "Npy2r", "Pvalb", "Nrn1"), nmols = 10000)
接下來，我們使用 SCTransform 進(jìn)行歸一化，然后進(jìn)行標(biāo)準(zhǔn)降維和聚類。此步驟從開始到結(jié)束大約需要 5 分鐘。
xenium.obj <- SCTransform(xenium.obj, assay = "Xenium")
xenium.obj <- RunPCA(xenium.obj, npcs = 30, features = rownames(xenium.obj))
xenium.obj <- RunUMAP(xenium.obj, dims = 1:30)
xenium.obj <- FindNeighbors(xenium.obj, reduction = "pca", dims = 1:30)
xenium.obj <- FindClusters(xenium.obj, resolution = 0.3)
然后，我們可以使用 DimPlot() 在 UMAP 空間中根據(jù)每個細(xì)胞的聚類對每個單元格進(jìn)行著色，或者使用 ImageDimPlot() 覆蓋在圖像上，從而可視化聚類結(jié)果。
DimPlot(xenium.obj)
我們可以在 UMAP 坐標(biāo)上可視化之前查看的標(biāo)記的表達(dá)水平。
FeaturePlot(xenium.obj, features = c("Cux2", "Bcl11b", "Foxp2", "Gad1", "Sst", "Gfap"))
現(xiàn)在，我們可以使用 ImageDimPlot() 為上一步中確定的簇標(biāo)簽著色的細(xì)胞位置著色。
ImageDimPlot(xenium.obj, cols = "polychrome", size = 0.75)
利用每個細(xì)胞的定位信息，我們能夠計算出它們各自的空間生態(tài)位。為了對細(xì)胞進(jìn)行注釋，我們采用了Allen Brain Institute提供的大腦皮層參考數(shù)據(jù)，因此首先需要將整個數(shù)據(jù)集的范圍限定在大腦皮層?？梢酝ㄟ^提供的鏈接安裝Allen Brain的參考數(shù)據(jù)。
在下面的分析中，我們利用Slc17a7基因的表達(dá)情況來輔助識別大腦皮層的具體區(qū)域。
xenium.obj <- LoadXenium("/brahms/hartmana/vignette_data/xenium_tiny_subset")
p1 <- ImageFeaturePlot(xenium.obj, features = "Slc17a7", axes = TRUE, max.cutoff = "q90")
p1
crop <- Crop(xenium.obj[["fov"]], x = c(600, 2100), y = c(900, 4700))
xenium.obj[["crop"]] <- crop
p2 <- ImageFeaturePlot(xenium.obj, fov = "crop", features = "Slc17a7", size = 1, axes = TRUE, max.cutoff = "q90")
p2
FindTransferAnchors函數(shù)能夠用來整合空間轉(zhuǎn)錄組數(shù)據(jù)集中的單個斑點數(shù)據(jù)。而Seurat v5版本進(jìn)一步支持了一種名為Robust Cell Type Decomposition（穩(wěn)健細(xì)胞類型分解，簡稱RCTD）的計算方法，該方法能夠在提供單細(xì)胞RNA測序（scRNA-seq）參考數(shù)據(jù)的基礎(chǔ)上，對空間數(shù)據(jù)集中的斑點數(shù)據(jù)進(jìn)行解卷積分析。RCTD已被證實能夠有效地對來自SLIDE-seq、Visium以及10x Genomics公司的Xenium原位空間平臺等多種技術(shù)的空間數(shù)據(jù)進(jìn)行注釋。
為了執(zhí)行RCTD分析，我們首先需要從GitHub上安裝名為spacexr的R包，該包提供了RCTD功能的實現(xiàn)。
devtools::install_github("dmcable/spacexr", build_vignettes = FALSE)
從 Seurat 查詢和參考對象中提取計數(shù)、聚類和點信息，以構(gòu)建 RCTD 用于注釋的參考和 SpatialRNA 對象。然后將注釋的輸出添加到 Seurat 對象。
library(spacexr)
query.counts <- GetAssayData(xenium.obj, assay = "Xenium", slot = "counts")[, Cells(xenium.obj[["crop"]])]
coords <- GetTissueCoordinates(xenium.obj[["crop"]], which = "centroids")
rownames(coords) <- coords$cell
coords$cell <- NULL
query <- SpatialRNA(coords, query.counts, colSums(query.counts))
# allen.corted.ref can be downloaded here:
# https://www.dropbox.com/s/cuowvm4vrf65pvq/allen_cortex.rds?dl=1
allen.cortex.ref <- readRDS("/brahms/shared/vignette-data/allen_cortex.rds")
allen.cortex.ref <- UpdateSeuratObject(allen.cortex.ref)
Idents(allen.cortex.ref) <- "subclass"
# remove CR cells because there aren't enough of them for annotation
allen.cortex.ref <- subset(allen.cortex.ref, subset = subclass != "CR")
counts <- GetAssayData(allen.cortex.ref, assay = "RNA", slot = "counts")
cluster <- as.factor(allen.cortex.ref$subclass)
names(cluster) <- colnames(allen.cortex.ref)
nUMI <- allen.cortex.ref$nCount_RNA
names(nUMI) <- colnames(allen.cortex.ref)
nUMI <- colSums(counts)
levels(cluster) <- gsub("/", "-", levels(cluster))
reference <- Reference(counts, cluster, nUMI)
# run RCTD with many cores
RCTD <- create.RCTD(query, reference, max_cores = 8)
RCTD <- run.RCTD(RCTD, doublet_mode = "doublet")
annotations.df <- RCTD@results$results_df
annotations <- annotations.df$first_type
names(annotations) <- rownames(annotations.df)
xenium.obj$predicted.celltype <- annotations
keep.cells <- Cells(xenium.obj)[!is.na(xenium.obj$predicted.celltype)]
xenium.obj <- subset(xenium.obj, cells = keep.cells)
與以往的獨立細(xì)胞分析不同，空間數(shù)據(jù)讓我們能夠更全面地理解細(xì)胞，不僅考慮它們周圍的局部環(huán)境，還包括它們在整個空間中的背景。在Seurat v5版本中，我們新增了對空間數(shù)據(jù)進(jìn)行“生態(tài)位”分析的功能，這種分析可以識別出組織中的特定區(qū)域（即“生態(tài)位”），每個區(qū)域都有其獨特的空間鄰近細(xì)胞類型組合。這一方法受到了Goltsev等人在2018年發(fā)表于《細(xì)胞》雜志和He等人在2022年發(fā)表于《自然生物技術(shù)》雜志的研究方法的啟發(fā)，我們?yōu)槊總€細(xì)胞定義了一個“局部鄰域”，這包括了與其在空間上距離最近的k個鄰居，并統(tǒng)計這個鄰域內(nèi)每種細(xì)胞類型的頻率。接著，我們應(yīng)用k均值聚類算法，將具有相似鄰域特征的細(xì)胞歸類到相同的空間生態(tài)位中。
在Seurat中，我們通過調(diào)用BuildNicheAssay函數(shù)來創(chuàng)建一個新的分析模塊，稱為“生態(tài)位”，它包含了每個細(xì)胞周圍空間的細(xì)胞類型組成信息。此外，該函數(shù)還返回了一個元數(shù)據(jù)列“niches”，該列展示了基于生態(tài)位分析的細(xì)胞聚類結(jié)果。
xenium.obj <- BuildNicheAssay(object = xenium.obj, fov = "crop", group.by = "predicted.celltype",
niches.k = 5, neighbors.k = 30)
然后，我們可以根據(jù)細(xì)胞類型身份或利基身份對細(xì)胞進(jìn)行分組。確定的生態(tài)位清楚地劃分了皮質(zhì)中的神經(jīng)元層。
celltype.plot <- ImageDimPlot(xenium.obj, group.by = "predicted.celltype", size = 1.5, cols = "polychrome",
dark.background = F) + ggtitle("Cell type")
niche.plot <- ImageDimPlot(xenium.obj, group.by = "niches", size = 1.5, dark.background = F) + ggtitle("Niches") +
scale_fill_manual(values = c("#442288", "#6CA2EA", "#B5D33D", "#FED23F", "#EB7D5B"))
celltype.plot | niche.plot
此外，我們觀察到每個生態(tài)位的組成對于不同的細(xì)胞類型都是豐富的。
table(xenium.obj$predicted.celltype, xenium.obj$niches)
##
## 1 2 3 4 5
## Astro 115 484 241 146 89
## Endo 27 250 66 62 45
## L2-3 IT 0 1749 14 5 6
## L4 2 2163 3 94 14
## L5 IT 2 66 2 627 84
## L5 PT 2 92 4 711 21
## L6 CT 82 2 0 34 857
## L6 IT 4 22 1 56 275
## L6b 95 0 0 2 2
## Lamp5 4 77 61 8 7
## Macrophage 7 48 15 16 10
## Meis2 0 0 0 0 0
## NP 0 1 0 78 9
## Oligo 305 130 42 76 70
## Peri 6 40 4 10 12
## Pvalb 5 155 0 75 54
## Serpinf1 0 5 0 1 1
## SMC 2 34 2 12 2
## Sncg 3 15 1 0 5
## Sst 0 53 0 55 26
## Vip 3 85 5 14 4
## VLMC 2 31 257 7 2
本文由mdnice多平臺發(fā)布
