对小鼠脑组织进行空间转录组学分析

生信
5.1k 词

和scRNA-seq比起来，就是每个样本点有了空间信息，并且每个样本点是一群细胞，那么要做的和之前的十分相似，只是主要有以下不同：

  • 寻找差异表达基因时，有2种手段
    • 聚类后在簇之间寻找，不需要空间信息
    • 根据空间信息寻找
  • 进行类型注释时，有2种手段
    • 寻找和scRNA-seq之间的锚点，然后整合，每个spot只能预测出1种细胞类型
    • 使用scRNA-seq作为参考，对每个spot进行反卷积，能得到每个spot内不同类型的比例

在这个示例里，都采用前者，原始代码见seurat官网。

数据导入

这里使用feature-barcode矩阵和空间图片作为输入。

数据导入 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
pkgs <- c("Seurat", "dplyr", "ggplot2", "patchwork")
pkgs_bioconductor <- c("glmGamPoi")
for(pkg in pkgs){
  if(!requireNamespace(pkg)) {
    install.packages(pkg, dependencies = TRUE)
  }
  library(pkg, character.only = TRUE)
}
for(pkg in pkgs_bioconductor){
  if(!requireNamespace(pkg)) {
    BiocManager::install(pkg, update = FALSE, ask = FALSE)
  }
  library(pkg, character.only = TRUE)
}

mtx_sp <- Read10X("./sourceData/filtered_feature_bc_matrix/")

image_sp <- Read10X_Image(image.dir = "./sourceData/spatial/",image.name = "tissue_hires_image.png")

# 详见https://github.com/satijalab/seurat/discussions/4833
image_sp@scale.factors$lowres = image_sp@scale.factors$hires

brain <- CreateSeuratObject(counts = mtx_sp, assay = "Spatial")

image_sp <- image_sp[colnames(brain)]

brain[['image']] <- image_sp

rm(list = setdiff(ls(), "brain"))
gc()

降维、聚类等

对于空转数据，相较于NormalizeData()的默认方法"LogNormalize"，用SCTransform效果更好。因为根据其步骤，LogNormalize会掩盖掉每个样本点总表达量的区别，而在空转数据中这很可能是有生物学意义的，比如这里的白质和肿瘤中坏死的区域可能呈现出这种特征。

在簇之间寻找差异表达基因，这里选择和原文一样的2个簇。

降维、聚类等 
1
2
3
4
5
6
7
8
9
10
11
12
brain <- SCTransform(brain, assay = "Spatial", verbose = FALSE)
brain <- RunPCA(brain, assay = "SCT", verbose = FALSE)

pcs_use <- (1:length(brain[["pca"]]@stdev))[cumsum(brain[["pca"]]@stdev^2) / sum(brain[["pca"]]@stdev^2) < 0.95]

brain <- FindNeighbors(brain, reduction = "pca", dims = pcs_use)
brain <- FindClusters(brain, verbose = FALSE)
brain <- RunUMAP(brain, reduction = "pca", dims = pcs_use)


de_markers <- FindMarkers(brain, ident.1 = 5, ident.2 = 7)
SpatialFeaturePlot(object = brain, features = rownames(de_markers)[1:3], alpha = c(0.1, 1), ncol = 3)

类型注释

这里先把空转数据切割成和教程差不多的形状，并在其中重新进行降维。

1
2
3
4
5
6
7
8
9
10
11
12
13
coords <- GetTissueCoordinates(brain)
brain$coords_x <- coords$x
brain$coords_y <- coords$y

cortex <- subset(brain, idents = c(4, 6, 3, 2,1,7))

SpatialFeaturePlot(cortex, features = c("coords_x","coords_y"))

cortex <- subset(cortex,subset = (coords_y < 6000 | coords_y < 19000 - 2*coords_x) & (coords_x < 7800 & coords_y > 3000))


cortex <- SCTransform(cortex, assay = "Spatial", verbose = FALSE) %>%
  RunPCA(verbose = FALSE)

然后和scRNA-seq数据整合。

1
2
3
4
5
6
7
8
9
10
11
12
13
allen_reference <- SCTransform(allen_reference, ncells = 3000, verbose = FALSE)

save(allen_reference,file = "allen_reference_SCTed.rdata")



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"

然后我们可以找出哪些细胞类型在空间上具有显著的变化。

SpatiallyVariableFeatures()好像有点问题

这是这个issue改进的代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
SpatiallyVariableFeatures_workaround <- function(object, assay="SCT", selection.method = "moransi") {
  #' This is work around function to replace SeuratObject::SpatiallyVariableFeatures function.
  #' return ranked list of Spatially Variable Features
  
  # Check if object is a Seurat object
  if (!inherits(object, "Seurat")) {
    stop("object must be a Seurat object")
  }
  
  # Check if assay is a valid assay
  if (!assay %in% names(object@assays)) {
    stop("assay must be a valid assay")
  }
  
  # Extract meta.features from the specified object and assay
  data <- object@assays[[assay]]@meta.features
  
  # Select columns starting with the provided col_prefix
  moransi_cols <- grep(paste0("^", selection.method), colnames(data), value = TRUE)
  
  # Filter rows where "moransi.spatially.variable" is TRUE
  filtered_data <- data[data[[paste0(selection.method, ".spatially.variable")]], moransi_cols]
  
  # Sort filtered data by "moransi.spatially.variable.rank" column in ascending order
  sorted_data <- filtered_data[order(filtered_data[[paste0(selection.method, ".spatially.variable.rank")]]), ]
  
  # Return row names of the sorted data frame
  rownames(sorted_data)
}
1
2
3
4
5
6
7
8
cortex <- FindSpatiallyVariableFeatures(cortex, assay = "predictions", selection.method = "moransi",
                                        features = rownames(cortex), r.metric = 5, slot = "data")
top.clusters <- head(SpatiallyVariableFeatures_workaround(cortex,assay = "predictions", selection.method = "moransi"), 6)


SpatialPlot(object = cortex, features = top.clusters, ncol = 3) & theme(legend.position = "right")
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)) & theme(legend.position = "right")

变化最显著的：

全部的：

额差不多吧。