作者，Evil Genius

人生多艰，且行且珍惜

这一篇我们来分享一个关于单细胞空间数据分析共表达网络的方法----hdWGCNA，文章在hdWGCNA identifies co-expression networks in high-dimensional transcriptomics data，2023年6月12日发表于Cell的子刊cell Reports methods,目前IF 8.109。

我们来解读一下并分享代码

生物系统是非常复杂的，在不同分子、细胞、器官和有机体之间严格调节的相互作用的基础上，被组织成一个多尺度的功能单元层次。虽然实验方法能够在数百万个细胞中进行转录组范围的测量，但流行的生物信息学工具不支持系统级分析。

hdWGCNA提供网络推理、基因模块识别、基因富集分析、统计测试和数据可视化等功能。hdWGCNA，一个单细胞和空间转录组学数据的共表达网络分析框架。共表达网络基于输入特征的转换成对相关性，从而对基因之间的相关性进行定量测量（原理和WGCNA类似）。

hdWGCNA的分析步骤

computing pairwise correlations of input features, weighting correlations with a
soft-power threshold（β）, computing the topological overlap between features, and unsupervised clustering via the Dynamic Tree Cut algorithm。单细胞数据中固有的稀疏性和噪声可能导致虚假的基因相关性，从而使共表达网络分析复杂化。此外，单细胞或空间转录组数据的相关结构在不同的亚群(细胞类型、细胞状态、解剖区域)中差异很大。scRNA-seq数据中的典型hdWGCNA工作流程通过将高度相似的细胞折叠成“metacell”来减少稀疏性，同时保留细胞异质性，并允许模块化设计在特定细胞群中执行单独的网络分析，从而考虑到这些因素。

metacell被定义为代表不同细胞状态的转录组相似细胞的小群。有几种方法可以从单细胞数据中识别metacell。hdWGCNA利用自引导聚合(bagging)算法，通过将k近邻(KNN)应用于输入数据集的降维表示，从单细胞数据集构建metacell转录组谱.

hdWGCNA从单细胞数据集和metacell表达矩阵中计算归一化基因表达矩阵中的基因-基因相关性，同时改变集成成单个metacell的细胞数量(KNN K参数)。在单细胞表达矩阵中，这些基因-基因相关性的分布在零处出现峰值，而在metacell矩阵中，扁平分布对应于更多的非零相关性，这表明与单细胞矩阵相比，metacell表达谱更不容易出现嘈杂的基因-基因相关性（降低噪音），单细胞本身数据的稀疏性也很大程度上降低。given a sufficient sample size, pseudo-bulk expression profiles are likely suitable for co-expression network analysis.

虽然共表达模块由许多基因组成，但将整个模块的表达汇总为一个指标是方便的。模块特征基因(MEs)被定义为模块基因表达矩阵的第一主成分，描述了整个共表达模块的表达模式。hdWGCNA使用高维数据的特定调整来计算MEs，允许批量校正和连续协变量的回归。另外，hdWGCNA还可以使用其他基因评分方法，如UCell或Seurat的AddModuleScore函数，分析发现这些评分与MEs相关。

然而，一些共表达模块可能对应于多种细胞类型共同的细胞过程，在这种情况下，hub基因可能被广泛表达。分析检查了选定的细胞类型特异性模块的MEs，发现其整体表达模式与其组成hub基因相似

使用ODC共表达网络展示了hdWGCNA的一些下游功能。为了实现网络可视化，使用统一流形近似和投影(UMAP)将共表达网络拓扑重叠矩阵(TOM)嵌入到二维流形中，使用每个基因的拓扑重叠与每个模块的top hub基因作为输入特征。发现10个ODC模块中有8个基于其MEs在ODC细胞中特异性表达。最后，进行了模块保存分析，以测试这些模块在独立数据集中的可重复性，并发现所有odc特定的模块都被显著保存。总之，这些人类PFC数据集中的网络分析显示了hdWGCNA工作流程的核心功能。

Spatial co-expression networks represent regional expression patterns

Spatial co-expression networks represent regional expression patterns in the mouse brain

细胞发育中的hdWGCNA分析

Isoform co-expression network analysis reveals fate-specific expression programs in the hippocampal radial glia lineage

方法注意：

单细胞数据形成metacell

空间数据形成metacell

Aggregation of neighboring spatial transcriptomic spots to form metaspots。

This approach leverages spatial coordinates rather than the dimensionality-reduced representation. For each ST spot, we obtain a list of physically neighboring spots. We then devise a grid of ’’principal spots’’, which are evenly spaced spots throughout the input tissue which serve as anchor points for aggregating neighboring spots. Each principal spot and its neighbors are aggregated into one metaspot, with at most seven spots merging into one metaspot and at most two overlapping spots between metaspots. We implemented this procedure as part of the hdWGCNA R package in mthe MetaspotsByGroups function. Similar to the MetacellsByGroups function, the user may specify groups within the Seurat object to perform the aggregation, such that metacells would only be grouped within the same tissue slice, anatomical region, or other annotation. For all downstream analysis with hdWGCNA, the metaspot expression dataset can be used in place of the metacell expression matrix.

计算共表达网络，与普通的WGCNA分析方法近乎一致。

最后来看看实例代码（high dimensional WGCNA）

• 安装

# create new conda environment for R
conda create -n hdWGCNA -c conda-forge r-base r-essentials

# activate conda environment
conda activate hdWGCNA

# install BiocManager
install.packages("BiocManager")

# install Bioconductor core packages
BiocManager::install()

# install additional packages:
install.packages(c("Seurat", "WGCNA", "igraph", "devtools"))

devtools::install_github('smorabit/hdWGCNA', ref='dev')

单细胞示例

# single-cell analysis package
library(Seurat)

# plotting and data science packages
library(tidyverse)
library(cowplot)
library(patchwork)

# co-expression network analysis packages:
library(WGCNA)
library(hdWGCNA)

library(igraph)

# using the cowplot theme for ggplot
theme_set(theme_cowplot())

# set random seed for reproducibility
set.seed(12345)

# optionally enable multithreading
enableWGCNAThreads(nThreads = 8)

# load the Zhou et al snRNA-seq dataset
seurat_obj <- readRDS('Zhou_2020.rds')

p <- DimPlot(seurat_obj, group.by='cell_type', label=TRUE) +
   umap_theme() + ggtitle('Zhou et al Control Cortex') + NoLegend()

p

seurat_obj <- SetupForWGCNA(
  seurat_obj,
  gene_select = "fraction", # the gene selection approach
  fraction = 0.05, # fraction of cells that a gene needs to be expressed in order to be included
  wgcna_name = "tutorial" # the name of the hdWGCNA experiment
)

Construct metacells

# construct metacells  in each group
seurat_obj <- MetacellsByGroups(
  seurat_obj = seurat_obj,
  group.by = c("cell_type", "Sample"), # specify the columns in seurat_obj@meta.data to group by
  reduction = 'harmony', # select the dimensionality reduction to perform KNN on
  k = 25, # nearest-neighbors parameter
  max_shared = 10, # maximum number of shared cells between two metacells
  ident.group = 'cell_type' # set the Idents of the metacell seurat object
)

# normalize metacell expression matrix:
seurat_obj <- NormalizeMetacells(seurat_obj)

seurat_obj <- SetDatExpr(
  seurat_obj,
  group_name = "INH", # the name of the group of interest in the group.by column
  group.by='cell_type', # the metadata column containing the cell type info. This same column should have also been used in MetacellsByGroups
  assay = 'RNA', # using RNA assay
  slot = 'data' # using normalized data
)

# Test different soft powers:
seurat_obj <- TestSoftPowers(
  seurat_obj,
  networkType = 'signed' # you can also use "unsigned" or "signed hybrid"
)

# plot the results:
plot_list <- PlotSoftPowers(seurat_obj)

# assemble with patchwork
wrap_plots(plot_list, ncol=2)

Construct co-expression network

# construct co-expression network:
seurat_obj <- ConstructNetwork(
  seurat_obj, soft_power=9,
  setDatExpr=FALSE,
  tom_name = 'INH' # name of the topoligical overlap matrix written to disk
)

PlotDendrogram(seurat_obj, main='INH hdWGCNA Dendrogram')

可视化

ModuleNetworkPlot(seurat_obj)

Combined hub gene network plots

# hubgene network
HubGeneNetworkPlot(
  seurat_obj,
  n_hubs = 3, n_other=5,
  edge_prop = 0.75,
  mods = 'all'
)

g <- HubGeneNetworkPlot(seurat_obj,  return_graph=TRUE)

Applying UMAP to co-expression networks

seurat_obj <- RunModuleUMAP(
  seurat_obj,
  n_hubs = 10, # number of hub genes to include for the UMAP embedding
  n_neighbors=15, # neighbors parameter for UMAP
  min_dist=0.1 # min distance between points in UMAP space
)

# get the hub gene UMAP table from the seurat object
umap_df <- GetModuleUMAP(seurat_obj)

# plot with ggplot
ggplot(umap_df, aes(x=UMAP1, y=UMAP2)) +
  geom_point(
   color=umap_df$color, # color each point by WGCNA module
   size=umap_df$kME*2 # size of each point based on intramodular connectivity
  ) +
  umap_theme()

ModuleUMAPPlot(
  seurat_obj,
  edge.alpha=0.25,
  sample_edges=TRUE,
  edge_prop=0.1, # proportion of edges to sample (20% here)
  label_hubs=2 ,# how many hub genes to plot per module?
  keep_grey_edges=FALSE

)

空间数据的示例

# 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 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)

    tissue row col imagerow imagecol
AAACAAGTATCTCCCA-1_1      1  50 102     7475     8501
AAACACCAATAACTGC-1_1      1  59  19     8553     2788
AAACAGAGCGACTCCT-1_1      1  14  94     3164     7950
AAACAGCTTTCAGAAG-1_1      1  43   9     6637     2099
AAACAGGGTCTATATT-1_1      1  47  13     7116     2375
AAACATGGTGAGAGGA-1_1      1  62   0     8913     1480

# 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

# add annotations to Seurat object
annotations <- read.csv('annotations.csv')
ix <- match(seurat_obj$seurat_clusters, annotations$seurat_clusters)
seurat_obj$annotation <- annotations$annotation[ix]

# set idents
Idents(seurat_obj) <- seurat_obj$annotation

p3 <- SpatialDimPlot(seurat_obj, label = TRUE, label.size = 3)
p3 + NoLegend()

构建metacell矩阵

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)

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')

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])

可视化

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

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()

当然，还有起他的分析内容

Network visualization

生活很好，有你更好