空间转录组第七讲：多样本合并、marker基因分析

前面介绍了利用Seurat对10x空间转录组数据单个样本进行数据预处理、降维、聚类分析。今天来继续分享怎么整合多个空间转录组样本，同时进行聚类和分析marker基因。跟单细胞一样，我们做空间转录组一般不太可能只做一个样本，一般会多个样本同时分析，样本之前往往需要进行亚群数目和基因表达差异的比较，这时候就需要把多个样本整合起来一起分析。Seurat同样提供空间转录组样本合并的函数merge，其实跟单细胞数据合并里的merge函数是一样的，就是不做批次效应去除直接整合在一起。Seurat推荐空间转录组数据使用sctransform 归一化方法，在一定程度上能减少样本间的批次效应的影响，但是如果组织切片的差异特别大，甚至不是同一种组织，建议还是不要放到一起合并分析啦。
上一节我们已经介绍了单个样本怎么读取和归一化，今天直接同样本合并开始。这里直接导入SeuratDataR包里内置小鼠脑的数据集，SeuratData是一种使用R的内部包和数据管理系统以Seurat对象的形式分发数据集的机制。它代表了一种方便用户访问Seurat中使用的数据集的方法。导入SeuratData R包后，我们可以用InstallData函数来安装下载对于的数据集。

SeuratData R包安装：

devtools::install_github('satijalab/seurat-data')
library(SeuratData)
AvailableData()

选择某个数据集，这里下载小鼠脑的数据集，选择anterior1和posterior1两个样本：

InstallData("stxBrain")
mydata1 <- LoadData("stxBrain", type = "anterior1")
mydata2 <- LoadData("stxBrain", type = "posterior1")

注意，如果是自己的样本，请使用上一讲单个样本处理的代码分别分析每个样本然到归一化那一步，然后进行后面的合并。

样本合并

data.merge <- merge(mydata1, mydata2)
DefaultAssay(data.merge) <- "SCT"
VariableFeatures(data.merge) <- c(VariableFeatures(mydata1), VariableFeatures(mydata2))

降维、聚类

data.merge <- RunPCA(data.merge, verbose = FALSE)
data.merge <- FindNeighbors(data.merge, dims = 1:30)
data.merge <- FindClusters(data.merge, verbose = FALSE)
data.merge <- RunUMAP(data.merge, reduction = "pca",dims = 1:30)
data.merge <- RunTSNE(data.merge, reduction = "pca",dims = 1:30)

umap、tsne可视化

p1 <- DimPlot(data.merge, reduction = "umap", label = TRUE)
p2 <- DimPlot(data.merge, reduction = "tsne", label = TRUE)
plot_grid(p1, p2)

image.png

切片组织亚群分布

SpatialDimPlot(data.merge, label = TRUE, label.size = 3)

image.png

不显示label

SpatialDimPlot(data.merge, label = FALSE, label.size = 3)

image.png

大脑分区marker基因表达分布展示

皮层marker基因STX1A

SpatialFeaturePlot(data.merge, features = c("Stx1a"))

image.png

丘脑marker基因Prkcd

SpatialFeaturePlot(data.merge, features = c("Prkcd"))

image.png

海马marker基因HPCA

SpatialFeaturePlot(data.merge, features = c("Hpca"))

image.png

脉络丛marker基因Ttr

SpatialFeaturePlot(data.merge, features = c("Ttr"))

image.png

也可以查看亚群中这几个marker基因分的表达分布

FeaturePlot(data.merge, features = c('Hpca','Ttr',"Prkcd","Stx1a"), cols = c("grey", "red"),reduction = "tsne")

image.png

亚群marker基因分析

同时分析所有亚群的marker基因

data.markers <- FindAllMarkers(data.merge, only.pos = FALSE, min.pct = 0.25, logfc.threshold = 0.25,test.use = "wilcox")

参数解析：
logfc.threshold: logfc阈值，阈值以上的基因才输出，默认0.25
only.pos：是否只输出正的marker，也就是只输出在该亚群中高表达的marker，默认是FALSE
test.use：算法选择，indAllMarkers提供8种差异分析算法，默认是wilcox
min.pct: 最低表达比例，指基因在两组细胞中至少有一组中有表达的细胞的比例等于或超过该阈值

查看每个亚群top3差异基因

topgene<-data.markers %>% group_by(cluster) %>% top_n(n = 3, wt = avg_logFC)
topgene
# A tibble: 69 x 7
# Groups:   cluster [23]
       p_val avg_logFC pct.1 pct.2 p_val_adj cluster gene
       <dbl>     <dbl> <dbl> <dbl>     <dbl> <fct>   <chr>
 1 0.            2.79  1     0.6   0.        0       Penk
 2 0.            2.78  1     0.357 0.        0       Gpr88
 3 0.            2.72  1     0.648 0.        0       Ppp1r1b
 4 2.69e-135     0.653 1     0.953 5.08e-131 1       Sparc
 5 1.51e-132     0.808 0.918 0.498 2.85e-128 1       Agt
 6 9.82e- 55     0.637 0.418 0.174 1.85e- 50 1       Tcf7l2
 7 0.            1.34  0.927 0.184 0.        2       Myl4
 8 0.            1.24  0.89  0.117 0.        2       Trbc2
 9 1.12e-307     1.64  0.998 0.521 2.12e-303 2       Vxn
10 1.21e-213     1.76  1     0.942 2.28e-209 3       Plp1
# … with 59 more rows

也可以单独两个或多个亚群比较分析marker基因分析

markers <- FindMarkers(data.merge, ident.1 = 2,, ident.2 = 6)

绘制每个亚群top3差异基因热图

DoHeatmap(data.merge, features = topgene$gene,size = 2) + NoLegend()

image.png

绘制marker基因小提琴图

VlnPlot(data.merge,features = c('Penk','Sparc'), group.by = "seurat_clusters",pt.size = 0)

image.png

绘制marker基因tsne图

FeaturePlot(data.merge, features = c('Penk','Sparc'), cols = c("grey", "red"),reduction = "tsne")

image.png

今天先分享到这，下期继续！

image.png

扫码关注公众号，内容更精彩哦！