前言
Immugent在hdWGCNA系列推文第一篇：hdWGCNA：将单细胞和空间转录组的WGCNA分析变成现实中,介绍了hdWGCNA包的主要功能框架，并在上一期推文:hdWGCNA系列推文（一）：分析流程的搭建中给大家说明了如何安装hdWGCNA。紧接着，在第三期推文：hdWGCNA系列推文（二）：分析单细胞转录组数据中介绍了如何用hdWGCNA包分析单细胞测序数据。那么今天，Immugent继续介绍hdWGCNA包的其它功能。
在上一篇推文中，Immugent已经介绍过hdWGCNA包最大的亮点就是其能进行单细胞测序数据的WGCNA分析。事实上，hdWGCNA包的厉害之处不仅仅在于其能分析单细胞数据，而且还能分析空间转录组数据。一般来说，通过分析单细胞数据的基因模块，我们可以找到一些规律，或者在一定程度上反应了生物学事件，但是无法真正意义上认定就是某些细胞的真正相互作用。因此，通过寻找在空间上的共定位信号，更能准确的反应真实性的细胞相互作用。
废话不多说，下面开始展示如何使用hdWGCNA包分析空间转录组数据。
代码流程
Load required libraries
# single-cell analysis package
library(Seurat)
# package to install the mouse brain dataset
library(SeuratData)
# plotting and data science packages
library(tidyverse)
library(cowplot)
library(patchwork)
# co-expression network analysis packages:
library(WGCNA)
library(hdWGCNA)
# install this package, which allows us to compute distance between the spots
install.packages('proxy')
library(proxy)
# enable parallel processing for network analysis (optional)
enableWGCNAThreads(nThreads = 8)
# using the cowplot theme for ggplot
theme_set(theme_cowplot())
# set random seed for reproducibility
set.seed(12345))
Download and process the mouse brain dataset
# download the mouse brain ST dataset (stxBrain)
SeuratData::InstallData("stxBrain")
# load the anterior and posterior samples
brain <- LoadData("stxBrain", type = "anterior1")
brain$region <- 'anterior'
brain2 <- LoadData("stxBrain", type = "posterior1")
brain2$region <- 'posterior'
# merge into one seurat object
seurat_obj <- merge(brain, brain2)
seurat_obj$region <- factor(as.character(seurat_obj$region), levels=c('anterior', 'posterior'))
# save unprocessed object
saveRDS(seurat_obj, file='mouse_brain_ST_unprocessed.rds')
# make a dataframe containing the image coordinates for each sample
image_df <- do.call(rbind, lapply(names(seurat_obj@images), function(x){
seurat_obj@images[[x]]@coordinates
}))
# merge the image_df with the Seurat metadata
new_meta <- merge(seurat_obj@meta.data, image_df, by='row.names')
# fix the row ordering to match the original seurat object
rownames(new_meta) <- new_meta$Row.names
ix <- match(as.character(colnames(seurat_obj)), as.character(rownames(new_meta)))
new_meta <- new_meta[ix,]
# add the new metadata to the seurat object
seurat_obj@meta.data <- new_meta
head(image_df)
# normalization, feature selection, scaling, and PCA
seurat_obj <- seurat_obj %>%
NormalizeData() %>%
FindVariableFeatures() %>%
ScaleData() %>%
RunPCA()
# Louvain clustering and umap
seurat_obj <- FindNeighbors(seurat_obj, dims = 1:30)
seurat_obj <- FindClusters(seurat_obj,verbose = TRUE)
seurat_obj <- RunUMAP(seurat_obj, dims = 1:30)
# set factor level for anterior / posterior
seurat_mouse_vis$region <- factor(as.character(seurat_mouse_vis$region), levels=c('anterior', 'posterior'))
# show the UMAP
p1 <- DimPlot(seurat_obj, label=TRUE, reduction = "umap", group.by = "seurat_clusters") + NoLegend()
p1
p2 <- SpatialDimPlot(seurat_obj, label = TRUE, label.size = 3)
p2
image.png
image.png
Construct metaspots
seurat_obj <- SetupForWGCNA(
seurat_obj,
gene_select = "fraction",
fraction = 0.05,
wgcna_name = "vis"
)
seurat_obj <- MetaspotsByGroups(
seurat_obj,
group.by = c("region"),
ident.group = "region",
assay = 'Spatial',
slot = 'counts'
)
seurat_obj <- NormalizeMetacells(seurat_obj)
m_obj <- GetMetacellObject(seurat_obj)
m_obj
image.png
Co-expression network analysis
# set up the expression matrix, set group.by and group_name to NULL to include all spots
seurat_obj <- SetDatExpr(
seurat_obj,
group.by=NULL,
group_name = NULL
)
# test different soft power thresholds
seurat_obj <- TestSoftPowers(seurat_obj)
plot_list <- PlotSoftPowers(seurat_obj)
wrap_plots(plot_list, ncol=2)
# construct co-expression network:
seurat_obj <- ConstructNetwork(
seurat_obj,
tom_name='test',
overwrite_tom=TRUE
)
# plot the dendrogram
PlotDendrogram(seurat_obj, main='Spatial hdWGCNA dendrogram')
image.png
image.png
seurat_obj <- ModuleEigengenes(seurat_obj)
seurat_obj <- ModuleConnectivity(seurat_obj)
seurat_obj <- ResetModuleNames(
seurat_obj,
new_name = "SM"
)
modules <- GetModules(seurat_obj) %>% subset(module != 'grey')
head(modules[,1:3])
Data visualization
# get module eigengenes and gene-module assignment tables
MEs <- GetMEs(seurat_obj)
modules <- GetModules(seurat_obj)
mods <- levels(modules$module); mods <- mods[mods != 'grey']
# add the MEs to the seurat metadata so we can plot it with Seurat functions
seurat_obj@meta.data <- cbind(seurat_obj@meta.data, MEs)
# plot with Seurat's DotPlot function
p <- DotPlot(seurat_obj, features=mods, group.by = 'annotation', dot.min=0.1)
# flip the x/y axes, rotate the axis labels, and change color scheme:
p <- p +
coord_flip() +
RotatedAxis() +
scale_color_gradient2(high='red', mid='grey95', low='blue') +
xlab('') + ylab('')
p
image.png
p <- SpatialFeaturePlot(
seurat_obj,
features = mods,
alpha = c(0.1, 1),
ncol = 8
)
p
png("figures/MEs_featureplot.png", height=16, width=20, units='in', res=200)
p
dev.off()
image.png
说在最后
从最后的分析结果中我们可以看出，使用hdWGCNA包不仅可以算出空间转录组的特征基因模块，还可以展示这些模块在空间上的分布。如果有对应的单细胞转录组数据，还可以利用反卷积算法，算出这些特定模块中准确的细胞类型，这样就能清晰的表征相互作用的细胞类型。此外，与传统使用差异分析算出的基因模块不同，运用hdWGCNA包算出的基因模块是具有一定功能特征的，而且都是高质量的功能性基因，这样就大大有助于我们做下游功能机制的研究。
好啦，本期分享到这里就结束了，我们下期再会~~
