江山代有套路出

大家好，我是晨曦，上次的推文给大家介绍了单细胞图谱类文章，相信大家不管是看过那篇推文，还是看了我们挑圈联靠其它单细胞的相关推文，对于单细胞，不管是从流程还是从分析方式上都应该不陌生了吧？

那么接下来晨曦就通过代码的形式，向大家展示单细胞空间转录组的分析流程，空间转录组绝对属于新的测序技术的新宠，2021年仅有160篇文章发表！如果说单细胞的福利是在前两年，那么当下的福利就理所应到的在空间转录组身上，紧跟分析红利，让我们一起走进单细胞空间转录的世界吧~

▌加载包

library(Seurat)library(SeuratData)library(ggplot2) library(patchwork)library(dplyr)#注意这里Seurat版本需要≥3.2

加载示例数据

InstallData("stxBrain")#如果代码下载困难,也可以尝试Seurat官网教程中提供了网页下载链接brain <- LoadData("stxBrain", type = "anterior1")View(brain)

#如果想要详细了解示例数据,可以采取在R中获取帮助文档的形式即'?stxBrain'

#如果采用网页下载,需要把下载数据加载到Seurat对象中,这一步需要使用函数'Load10×_Spatial()' #空间数据存储在Seurat的形式

1. A spot by gene expression matrix(表达矩阵)

2. An image of the tissue slice (obtained from H&E staining during data acquisition)(组织切片图像)

3. Scaling factors that relate the original high resolution image to the lower resolution image used here for visualization.

晨曦解读

广义来说,空间转录组学有两种方法,其中一种是在测序的时候,以一种保留空间信息的方式捕获mRNA,这里说一个报道可能会更方便大家理解,在2018年,10×公司收购了Spatial Transcriptomics公司,而这个被收购的公司开发的空间转录组技术能够从一片完整的冰冻切片中,获取切片上不同位置细胞中的完整转录组数据的这项技术自然到了10×公司手中,随后10×公司把这项技术加以提升,使捕获单元(spot)内有1-10个细胞(的基因及其表达量),而未来展位则是达到彻底的精确,即一个捕获单元内只有一个细胞(这里的spot我们可以理解为是一个人为制造"细胞")(全称对比scRNA-seq会更方便理解)

▌数据分析过程(对比scRNA-seq流程发现后缀添加了_Spatial))

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)

#首先让我们看一下数据,对于空间数据集,分子计数的差异可能很大,特别是如果整个组织的细胞密度存在差异,我们在这里就看到了大量的异质性,这就需要进行有效的标准化

# 同时,下面的图表明,不同的计数差异不仅是技术性的,而且还取决于组织解剖结构,再者如果单纯使用LogNormalize()函数强制每个都具有相同基础来标准化,可能会出现违反生物学背景的问题,所以建议使用SCTransofrm()函数来进行标准化

#数据标准化brain <- SCTransform(brain, assay= "Spatial", verbose= FALSE)#这里做一下简单的探索(探索lognormalized和SCTransform标准化的效益)# rerun normalization to store sctransform residualsfor allgenesbrain <- SCTransform(brain, assay = "Spatial", return.only.var.genes = FALSE, verbose = FALSE)# also runstandard log normalization for comparisonbrain <- NormalizeData(brain, verbose = FALSE, assay = "Spatial")# Computes thecorrelation of the log normalized data and sctransform residuals with thenumber 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("LogNormalization")+

theme(plot.title = element_text(hjust = 0.5))p2 <- GroupCorrelationPlot(brain, assay = "SCT", cor = "nCount_Spatial_cor") + ggtitle("SCTransformNormalization") + theme(plot.title = element_text(hjust = 0.5)) p1 + p2

#官方在标准化之前并没有进行数据质控(QC),个人觉得应该加上,毕竟QC是所有数据分析的开始,QC的流程可以参考scRNA-seq,无外乎就是过滤掉少基因的spot和每个基因至少在n个spot上表达(标准可以根据实际情况进行更改)还有可以去除核糖体基因的影响

晨曦解读

每个特征基因与UMI数量的相关性,根据基因的平均表达进行分组,并生成相关性的箱线图,可以很清晰的看到LogNormalization的归一化效能不如SCTransform(即归一化的目的在于并不是平衡数量,而是就算每一个位置的细胞数量不同,但是表达量和UMI的增进趋势应该是大致相同的)(简单举例子就是小明有很多女性朋友,应该一样亲密才不容易翻车,但如果一个过于亲密就就很容易翻车)

基因表达可视化

SpatialFeaturePlot(brain, features = c("Hpca", "Ttr"))

#这里的Hcpa和Ttr是Gene,类似与scRNA-seq展示单基因在亚群的表达情况,只不过这里换成了空间图像

#可以通过参数来调节图像,快来个性化DIY吧~

p1 <- SpatialFeaturePlot(brain, features = "Ttr", pt.size.factor = 1)#pt.size.factor参数调节点大小，默认值为1.6

p2 <- SpatialFeaturePlot(brain, features = "Ttr", alpha = c(0.1, 1))#alpha参数调节点的最小和最大透明度，默认值为c（1，1），尝试设置c（0.1，1）以降低具有较低表达的点的透明度

p1 + p2

数据降维\聚类\可视化

#对表达数据走常规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()把聚类亚群叠加在空间上 p1 <- DimPlot(brain, reduction = "umap", label = TRUE)#这个标签参数挺有意思,个人觉得很可爱~ p2 <- SpatialDimPlot(brain, label = TRUE, label.size = 3)

p1 + p2

# 使用cells.highlight参数高亮感兴趣的一些细胞

SpatialDimPlot(brain, cells.highlight = CellsByIdentities(object = brain, idents = c(1, 2, 5, 3, 4, 8)), facet.highlight = TRUE, ncol = 3)

交互式可视化(其实就是spatialRNA的可视化)

SpatialDimPlot(brain, interactive = TRUE)#您可以在其中将鼠标悬停在点上并查看单元名称和当前标识类SpatialFeaturePlot(brain, features = "Ttr", interactive = TRUE)

#对于 SpatialFeaturePlot()，将interactive 设置为TRUE 会显示一个交互式窗格，您可以在其中调整点的透明度、点大小以及正在绘制的检测和特征

#探索数据后，选择完成按钮将返回最后一个活动图作为 ggplot 对象。LinkedDimPlot(brain)

#LinkedDimPlot（）函数将UMAP表示链接到组织图像表示，并允许交互式选择。例如，您可以在 UMAP 图中选择一个区域，图像表示中的相应点将被突出显示

#这里介绍几个参数

#cells.highlight参数:传递spot名称向量,可以在spatialdimplot上突出显示#facet.highlight参数:逻辑值,按分组分面显示spot图像

鉴定空间可变基因

Seurat提供的两种方法来识别marker基因,一种仍然是我们scRNA-seq的常规流程聚类→细胞类群之间差异找marker基因(这样我们只能通过算法以及表达量来确定亚群或者基因所处的空间位置),另一种同时因为我们的样本不再是细胞而是一个组织切片,可以通过不同的解剖区域(这个解剖区域可以通过对不同基因表达量的探索得出)进行注释后,找寻不同区域之间的marker基因

#组织内预先注释解剖区域进行差异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)#下图展示的Gene显示出来明显的空间限制

#没有预先注释鉴定出不同空间模式（spatial patterning）brain <- FindSpatiallyVariableFeatures(brain, assay = "SCT", features = VariableFeatures(brain)[1:1000],selection.method = "markvariogram")#可视化前6个表达

top.features <- head(SpatiallyVariableFeatures(brain, selection.method = "markvariogram"), 6) SpatialFeaturePlot(brain, features = top.features, ncol = 3, alpha = c(0.1, 1))

可视化区域的子集(其实就是对比scRNA-seq亚群再细分)

#提取子集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

联合scRNA-seq数据

scRNA-seq和spatialRNA-seq能否联合,这里说一下个人观点,应该划分研究目的,如果目的是构建大型的细胞图谱,也就是广谱类研究,那么不同来源\不同测序方法的样本整合在一起是有意义的,但是如果你研究的疾病存在很大的异质性,就比如说肿瘤,你A患者用spatialRNA-seq,B患者用scRNA-seq,就好像把水果派和沙拉一起放在餐桌上,除了可以证明这些都是吃的以外(患者数据),还能说明什么呢?

这里有联合数据有两种方法,分别是CCA以及MIA

#首先读取数据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 cortexcortex <- 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

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#数据可视化DefaultAssay(cortex) <- "predictions"SpatialFeaturePlot(cortex, features = c("L2/3 IT", "L4"), pt.size.factor = 1.6, ncol = 2, crop = TRUE)

#整体逻辑就是首先提取spatialRNA数据中我们感兴趣的解剖区域的数据，这里我们选择cortex（皮质）的spot子集，然后重新对这个子集进行标准化#其次对cortex这个区域对应的单细胞数据进行读取并标准化

#通过寻找锚点将scRNA-seq与spatialRNA数据联系起来，并把scRNA-seq数据的cluster定位在组织切片上#基于预测分数,我们还可以预测特定解剖位置的细胞类型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)

使用Seurat处理多个切片数据(也就是处理多个空间转录组数据)

brain2 <- LoadData("stxBrain", type = "posterior1")brain2 <- SCTransform(brain2, assay = "Spatial", verbose = FALSE)# 使用 merge 函数合并多个 slices 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"))