重温seurat官方教程四：空间转录组分析
上次系统学习seurat教程还是它升级到V3版本的时候，一段时间不见发现官网又发布了不少新的教程，很多老的教程也重新编辑过。单细胞测序技术的应用越来越普及，大家对它的理解也越发深入，分析方法和流程也有相应的调整。Seurat作为单细胞分析最普及的工具，其作者对单细胞数据分析的见解是我们望尘莫及的，因此重温大佬的教程或许能更新我们的认知。此专辑并不忠于原文，会有很多我自己的想法。
往期专题
单细胞转录组基础分析专题
10X空间转录组简介
技术原理
10x空间转录组（spatialRNA）的技术类似于表达谱芯片，不同之处在于表达谱芯片上是cDNA杂交探针，空间转录组芯片上则是mRNA的逆转录引物。这些引物和scRNA技术一样，具有特殊的barcode和UMI序列，他们在芯片上成簇“种植”，10x公司把这些“簇”称为spot。每个spot中的引物都具有唯一且相同的barcode序列，因此后续分析可以通过barcode序列确定spot在芯片上的位置。冰冻组织切片覆盖在空间转录组芯片上之后，使用特殊的试剂让RNA分子从细胞里释放出来，这些RNA到达spot之后就可以与其中的引物杂交，逆转录之后就可以得到带有barcode和UMI序列的cDNA，然后就可以建库测序了。
空间转录组可以给我们什么呢？一张HE染色的冰冻组织切片图像，覆盖这张切片的几千个spot的转录组数据。
spatialRNA的优势
同时提供组织切片的HE染色图像和转录组数据，让我们可以使用图像和分子相互印证分析结果。
把空间位置信息和转录组分析结合，展示组织中不同位置的基因活跃信息。
没有单细胞悬浮液制备的步骤，实验难度相比scRNA-seq大大降低。组织限制方面，除了皮肤、胰腺和骨头组织外，其余样本类型均可开展空间转录组研究。对于scRNA-seq难度很大的神经组织和肌肉组织，spatialRNA可以轻松完成实验。
空间位置重要性的示例
下图A和B是两个肿瘤组织的切片，A样本中存在淋巴细胞浸润现象，预后效果良好；而B样本淋巴细胞都被肿瘤组织阻碍在肿瘤的边缘，预后效果差。
如果只是通过常规转录组或者单细胞转录组，研究人员在两个样本中均可以发现淋巴细胞相关基因或者淋巴细胞群，无法解释两个样本预后之间的差异。但是通过空间转录组，就能很清晰看到两个样本之间的差异。
创建spatialRNA对象
我们使用10x官网提供的小鼠大脑切片数据。Seurat教程使用SeuratData函数导入数据，考虑到大家分析数据的实际情况，我从cellranger的输出结果导入数据。
注意：此教程要求seurat>=3.2版本
##==创建spatialRNA对象==####设置目录dir.create("Vignette04")setwd("Vignette04")dir.create("Anterior")dir.create("Posterior")##下载Mouse Brain Serial Section 1 (Sagittal-Anterior)数据#表达矩阵download.file(url="https://cf.10xgenomics.com/samples/spatial-exp/1.1.0/V1_Mouse_Brain_Sagittal_Anterior/V1_Mouse_Brain_Sagittal_Anterior_filtered_feature_bc_matrix.h5", destfile = "Anterior/matrix.h5")#染色图像download.file(url="https://cf.10xgenomics.com/samples/spatial-exp/1.1.0/V1_Mouse_Brain_Sagittal_Anterior/V1_Mouse_Brain_Sagittal_Anterior_spatial.tar.gz", destfile = "Anterior/image.tar.gz")untar("Anterior/image.tar.gz", exdir="Anterior")#cloupefiledownload.file(url="https://cf.10xgenomics.com/samples/spatial-exp/1.1.0/V1_Mouse_Brain_Sagittal_Anterior/V1_Mouse_Brain_Sagittal_Anterior_cloupe.cloupe", destfile = "Anterior/Anterior.cloupe")##下载Mouse Brain Serial Section 1 (Sagittal-Posterior)数据#表达矩阵download.file(url="https://cf.10xgenomics.com/samples/spatial-exp/1.1.0/V1_Mouse_Brain_Sagittal_Posterior/V1_Mouse_Brain_Sagittal_Posterior_filtered_feature_bc_matrix.h5", "Posterior/matrix.h5")#染色图像download.file(url="https://cf.10xgenomics.com/samples/spatial-exp/1.1.0/V1_Mouse_Brain_Sagittal_Posterior/V1_Mouse_Brain_Sagittal_Posterior_spatial.tar.gz", destfile = "Posterior/image.tar.gz")untar("Posterior/image.tar.gz", exdir="Posterior")#cloupefiledownload.file(url="https://cf.10xgenomics.com/samples/spatial-exp/1.1.0/V1_Mouse_Brain_Sagittal_Posterior/V1_Mouse_Brain_Sagittal_Posterior_cloupe.cloupe", destfile = "Posterior/Posterior.cloupe")##创建spatialRNA对象library(Seurat)library(tidyverse)library(patchwork)brain <- Load10X_Spatial(data.dir="Anterior", filename="matrix.h5", slice="Anterior")brain2 <- Load10X_Spatial(data.dir="Posterior", filename="matrix.h5", slice="Posterior")
Load10X_Spatial()函数创建Spatial Seurat对象比较方便，它需要cellranger输出文件中的HDF5格式表达矩阵和保存图像的spatial文件夹。data.dir是用来指定保存文件的目录，filename是用来指定表达矩阵的文件名。
数据预处理
官方教程没有去除零表达的基因，没有过滤基因数少的spot，也没有去除核糖体基因对分析的影响，我对此做了一定的调整。
##==数据预处理==##setwd('Anterior')##按每个spot最少200基因，每个基因最少有一个spot表达的标准过滤brain <- stRNA_filter(brain, min.features=1, min.spots=200) #stRNA_filter是内部函数，非公开资源brain[["percent.mt"]] <- PercentageFeatureSet(brain, pattern = "^mt-")p1 <-VlnPlot(brain, features = c("nFeature_Spatial", "nCount_Spatial", "percent.mt"),pt.size = 0.1, ncol = 3) + NoLegend()p2 <-VlnPlot(brain, features = "nFeature_Spatial", pt.size = 0.1) + NoLegend() +theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank())p3 <- SpatialFeaturePlot(brain, features = "nCount_Spatial") + theme(legend.position = "right")p <- p2|p3ggsave("qc_mesures.png", p1, width=8, height=5)ggsave("qc_spatial.png", p, width=9, height=5)##数据标准化并去除核糖体的影响brain <- SCTransform(brain, assay = "Spatial", vars.to.regress="percent.mt")
空间转录组检测到基因数明显比scRNA-seq多一些
降维聚类
brain <- RunPCA(brain, assay = "SCT", verbose = FALSE)plot1 <- DimPlot(brain, reduction = "pca", group.by="orig.ident", label=T) + NoLegend()plot2 <- ElbowPlot(brain, ndims=50, reduction="pca")plotc <- plot1+plot2ggsave("pca.png", plot = plotc, width = 8, height = 4)
nPCs=1:40brain <- FindNeighbors(brain, reduction = "pca", dims = nPCs)brain <- FindClusters(brain, verbose = FALSE)brain <- RunTSNE(brain, reduction = "pca", dims = nPCs)brain <- RunUMAP(brain, reduction = "pca", dims = nPCs)p1 <- DimPlot(brain, reduction = "umap", label = TRUE)p2 <- SpatialDimPlot(brain, label = TRUE, label.size = 3)p3 <- DimPlot(brain, reduction = "tsne", label = TRUE)plotc1 <- p1|p2plotc2 <- p1|p3ggsave("SpatialDimPlot.png", plotc1, width=8, height=4)ggsave("UMAP_tSNE.png", plotc2, width=8, height=4)
spatialRNA可视化
Seurat支持在组织切片上原位展示spot的转录特征，它通过SpatialDimPlot、SpatialFeaturePlot、ISpatialDimPlot、ISpatialFeaturePlot、LinkedDimPlot等函数实现可视化，后三个函数是基于shiny的互动可视化图形。这些函数有些参数需要注意：
crop：通用参数，逻辑值，是否修剪原始图像以便集中显示选定的spot
interactive：通用参数，逻辑值，是否打开shiny互动图形，当前版本有bug
pt.size.factor：通用参数，调整spot大小的参数，默认值为1.6
features：SpatialFeaturePlot参数，需要在原位spot图上展示的features
alpha：SpatialFeaturePlot参数，spot透明度区间，seurat会将features的值映射为透明度
group.by：SpatialDimPlot参数，分组依据的meta.data列名
cells.highlight：SpatialDimPlot参数，传递spot名称向量，可以在spatialdimplot上突出显示
facet.highlight：SpatialDimPlot参数，逻辑值，按分组分面显示spot图像
##==spatialRNA可视化==####SpatialFeaturePlotp1 <- SpatialFeaturePlot(brain, features = "Hpca")p2 <- SpatialFeaturePlot(brain, features = "Hpca", crop=F)p3 <- SpatialFeaturePlot(brain, features = "Hpca", pt.size.factor=1.2)p4 <- SpatialFeaturePlot(brain, features = "Hpca", alpha=c(0,1))p = (p1|p2)/(p3|p4)ggsave("spot_feature_comparison.png", p, width=8, height=8)#互动SpatialFeaturePlotSpatialFeaturePlot(brain, features = "Hpca", interactive=T)#有bug，报错: Cannot find 'NA' in this Seurat object#替代代码ISpatialFeaturePlot(brain, feature = "Hpca")
个人觉得图1最好看
##SpatialDimPlotp1 <- SpatialDimPlot(brain, label = TRUE, label.size = 3) + ggtitle("p1") + NoLegend()p2 <- SpatialDimPlot(brain, cells.highlight = colnames(subset(brain, idents = c(1, 2)))) + ggtitle("p2") + NoLegend()p3 <- SpatialDimPlot(brain, cells.highlight = CellsByIdentities(object = brain, idents = c(1, 2, 4, 5, 7, 8)), facet.highlight = TRUE, ncol = 3)p = (p1|p2)/p3ggsave("spot_dim_comparison.png", p, width=6, height=8)#互动SpatialDimPlotSpatialDimPlot(brain, interactive = TRUE)#有bug，报错: Cannot find 'NA' in this Seurat object#替代代码ISpatialDimPlot(brain)
#展示umap与spot图之间联系的shiny应用LinkedDimPlot(brain)这条命令会打开一个shiny应用，你在左图选择的spot，右图也会同步高亮显示
识别空间marker基因
在scRNA-seq的数据分析中，我们一般使用无监督聚类将细胞分成多个clusters，然后寻找clusters之间的marker基因，最后通过这些marker基因识别clusters的细胞类型。空间转录组中，我们依然可以使用这套分析流程，但是我们还有其他选择。别忘了我们的样本是一张组织切片，专业人士可以通过染色图像识别不同的解剖区域，他们甚至可以使用SpatialFeaturePlot提供的基因表达情况来印证自己的判断。因此，我们在执行差异分析之前就可以对各个解剖区域进行注释，然后再寻找可以代表各个区域的更多的biomarker。
#==空间marker基因==###基于cluster或先验知识注释de_markers <- FindMarkers(brain, ident.1 = 4, ident.2 = 6)p <- SpatialFeaturePlot(object = brain, features = rownames(de_markers)[1:3], ncol = 3)ggsave("spatial_biomarkers.png", p, width=8, height=4)
在没有预先注释的情况下寻找空间marker基因，可以使用算法识别表达水平取决于其空间位置的基因。
brain <- FindSpatiallyVariableFeatures(brain, assay = "SCT",features = VariableFeatures(brain), selection.method = "markvariogram")top.features <- head(SpatiallyVariableFeatures(brain, selection.method = "markvariogram"), 3)p <- SpatialFeaturePlot(brain, features = top.features, ncol = 3, alpha = c(0.1, 1))ggsave("spatial_markvariogram.png", p, width=8, height=4)
整合单细胞数据
经常有人问我空间转录组和单细胞转录组数据可以整合吗，整合的数据可以来自不同的样本吗？scRNA与spatialRNA的整合目前有两个比较权威的算法，一个是seurat提供的基于锚点的整合算法，另一个是发表在Nat Biotechnol上的MIA算法。至于整合的数据能否来自不同样本，要看你的研究材料和目的。如果你研究目的是构建正常组织的单细胞图谱，那么不同样本来源的数据整合在一起是有意义的。如果你研究的材料在个体之间存在明显的异质性，比如上文提到的淋巴浸润的肿瘤组织样本，你用A样本的spatialRNA数据和B样本的scRNA数据整合，准备说明什么问题呢？
本文将按seurat官网教程介绍锚点整合的方法，所用单细胞数据来自 Allen Institute的成年小鼠皮质细胞数据集，包含使用SMART-Seq2技术测序得到的大约14000个细胞的转录组数据。数据下载链接：https://www.dropbox.com/s/cuowvm4vrf65pvq/allen_cortex.rds?dl=1
##提取cortex区域的spot子集cortex <- subset(brain, idents = c(0, 3, 4, 5, 7))cortex <- subset(cortex, anterior_imagerow > 400 | anterior_imagecol < 150, invert = TRUE)cortex <- subset(cortex, anterior_imagerow > 275 & anterior_imagecol > 370, invert = TRUE)cortex <- subset(cortex, anterior_imagerow > 250 & anterior_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)p <- p1 + p2ggsave("spatial_cortex.png", p, width=8, height=4)#按官方教程对cortex子集重新标准化cortex <- SCTransform(cortex, assay = "Spatial", verbose = FALSE) %>% RunPCA(verbose = FALSE)##读取Allen Institute的cortex区域单细胞数据集allen_cortex <- readRDS("allen_cortex.rds")#note that setting ncells=3000 normalizes the full dataset but learns noise models on 3k cells，只为运算加速allen_cortex <- SCTransform(allen_cortex, ncells = 3000, verbose = FALSE) %>% RunPCA(verbose = FALSE) %>% RunUMAP(dims = 1:30)anchors <- FindTransferAnchors(reference = allen_cortex, query = cortex, normalization.method = "SCT")predictions.assay <- TransferData(anchorset = anchors, refdata = allen_cortex$subclass,prediction.assay = TRUE, weight.reduction = cortex[["pca"]])cortex[["predictions"]] <- predictions.assaysubclass <- TransferData(anchorset = anchors, refdata = allen_cortex$subclass, weight.reduction = cortex[["pca"]])cortex <- AddMetaData(cortex, subclass[,1], col.name = 'subclass')p1 <- DimPlot(allen_cortex_sct, group.by = 'subclass', reduction = "umap", label = T) + NoLegend()p2 <- SpatialDimPlot(cortex, group.by = 'subclass', label = T, label.size = 2, repel = T) + NoLegend()p = p1|p2ggsave("integrated_dimplot.png", p, width=8, height=4)
通过锚点我们将scRNA与spatialRNA数据联系起来了，并把scRNA数据的clusters定位在组织切片上了。
DefaultAssay(cortex) <- "predictions"p <- SpatialFeaturePlot(cortex, features = c("Astro", "L2/3 IT", "L4", "L5 PT", "L5 IT", "L6 CT", "L6 IT", "L6b", "Oligo"), pt.size.factor = 1, ncol = 3, crop = FALSE, alpha = c(0.1, 1))ggsave("integrated_features.png", p, width=6, height=6)
我们可以把每个subclass的预测分值展现在组织切片上，并通过调整透明度参数alpha，让得分低的spot能见度降低。
合并多个spatialRNA数据
我们此前下载了两个空间转录组的数据，除了brain对象我们还创建了brain2对象，我们现在把两个切片的数据合并在一起。
brain2 <- SCTransform(brain2, assay = "Spatial", verbose = FALSE)brain.merge <- merge(brain, brain2)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)DimPlot(brain.merge, reduction = "umap", group.by = c("ident", "orig.ident"))SpatialDimPlot(brain.merge)p <- SpatialFeaturePlot(brain.merge, features = c("Hpca", "Plp1"))ggsave("spatial_merge.png", p, width=6, height=6)
如果您阅读此文有所疑惑，或有不同见解，亦或其他问题，欢迎添加我的微信探讨。如果您有相关的测序或分析项目，也欢迎与我洽谈合作。
优惠活动
我服务的公司与10x Genomics公司共同举办了10x Genomics中国区Visium空间科学挑战赛。此项活动奖励丰厚，一等奖价值20万元，有兴趣的朋友可以联系我报名参赛。