那么接下来晨曦就通过代码的形式,向大家展示单细胞空间转录组的分析流程,空间转录组绝对属于新的测序技术的新宠,2021年仅有160篇文章发表!如果说单细胞的福利是在前两年,那么当下的福利就理所应到的在空间转录组身上,紧跟分析红利,让我们一起走进单细胞空间转录的世界吧~ ▌加载包 library(Seurat) library(SeuratData) library(ggplot2) library(patchwork) library(dplyr) #注意这里Seurat版本需要≥3.2 加载示例数据
#如果想要详细了解示例数据,可以采取在R中获取帮助 文档的形式即'?stxBrain' #如果采用网页下载,需要把下载数据加载到Seurat对象中,这一步需要使用函数'Load10×_Spatial()' #空间数据存储在Seurat的形式 晨曦解读 ▌数据分析过程(对比scRNA-seq流程发现后缀添加了_Spatial))
#数据标准化 brain <- SCTransform(brain, assay= 'Spatial', verbose= FALSE) #这里做一下简单的探索(探索lognormalized和SCTransform标准化的效益) # rerun normalization to store sctransform residualsfor allgenes brain <- SCTransform(brain, assay = 'Spatial', return.only.var.genes = FALSE, verbose = FALSE) # also runstandard log normalization for comparison brain <- 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的增进趋势应该是大致相同的)(简单举例子就是小明有很多女性朋友,应该一样亲密才不容易翻车,但如果一个过于亲密就就很容易翻车) 基因表达可视化
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 数据降维\聚类\可视化
# 使用cells.highlight参数高亮感兴趣的一些细胞 SpatialDimPlot(brain, cells.highlight = CellsByIdentities(object = brain, idents = c(1, 2, 5, 3, 4, 8)), facet.highlight = TRUE, ncol = 3)
鉴定空间可变基因 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显示出来明显的空间限制
可视化区域的子集(其实就是对比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数据 这里有联合数据有两种方法,分别是CCA以及MIA
#识别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)
使用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'))
SpatialFeaturePlot(brain.merge, features = c('Hpca', 'Plp1')) 欢迎大家关注解螺旋生信频道-挑圈联靠公号~ |
|