0. 数据准备
- 输入数据集的要求：已经进行了如下分析的Seurat对象
- 导入演示数据
#官方演示数据集
wget https://swaruplab.bio.uci.edu/public_data/Zhou_2020.rds
seurat_obj <- readRDS('Zhou_2020.rds')
这是一个正常的脑组织数据集，包含了使用Harmony整合的12个样本的Seurat对象。
- 加载R包和数据集
# 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)
# using the cowplot theme for ggplot
theme_set(theme_cowplot())
# set random seed for reproducibility
set.seed(12345)
#可视化一下看看
p <- DimPlot(seurat_obj, group.by='cell_type', label=TRUE) +
umap_theme()
p
1. Set up Seurat object for WGCNA
在进行hdWGCNA运算之前，我们首先需要处理Seurat object。绝大多数 hdWGCNA运算所需要的信息都储存在Seurat object’s @misc slot。一个Seurat object can hold multiple hdWGCNA experiments, for example representing different cell types in the same single-cell dataset. Notably, since we consider hdWGCNA to be a downstream data analysis step, we do not support subsetting the Seurat object after SetupForWGCNA has been run.
在这一步我们使用 SetupForWGCNA 来指定 hdWGNCA experiment的名称。这一步指定了用于WGCNA分析的基因 。关于基因的选择，有如下三种方法：
在这个演示里，我们选择了至少在5%的细胞中表达的基因（3857个），并将hdWGCNA experiment 命名为 “tutorial”。
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
)
length(seurat_obj@misc$tutorial$wgcna_genes)
# [1] 12639
2. Construct metacells
准备好 Seurat 对象后，首先需要从单细胞数据集中构建metacells。 简而言之，metacells是来自相同生物样品的一小组相似细胞的集合。 k-最近邻 (KNN) 算法可以识别相似细胞组，然后计算这些细胞的平均或总表达量，从而产生metacells基因表达矩阵。 metacells表达矩阵的稀疏性与原始表达矩阵相比大大降低，因此更适用于后续分析（WGCNA 等相关网络方法对数据稀疏性很敏感）。（单细胞表观基因组方法，例如 Cicero，在构建可访问网络之前采用了类似的元细胞聚合方法。）
hdWGCNA 使用 MetacellsByGroups函数来构造给定单细胞数据集的metacells表达矩阵。 此函数为存储在 hdWGCNA experiment中的metacells数据集构造一个新的 Seurat 对象。 group.by 参数确定将在哪些组中构建metacells。我们只想从来自相同生物样本的细胞构建metacells，因此通过 group.by 参数将该信息传递给 hdWGCNA 至关重要。 此外，我们通常为每种细胞类型分别构建metacells。 因此，在这个演示中，我们按 Sample 和 cell_type 进行分组以达到预期的结果。
k值（The number of cells to be aggregated）应根据输入数据集的大小进行调整，通常较小的 k 数可用于小数据集。 我们一般使用 20 到 75 之间的 k 值。本教程使用的数据集有 40,039 个细胞，每个生物样本中的细胞数从 890 到 8,188 不等，这里我们使用 k=25。 可以使用 max_shared 参数调整元单元之间允许的重叠量。
注意：metacells聚合方法对extremely underrepresented cell types效果不好。 例如，在这个数据集中，脑血管细胞（周细胞和内皮细胞）是the least represented，因此将它们排除在此分析之外。 MetacellsByGroups 有一个参数 min_cells 来排除小于指定单元数的组。
# 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
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)
Optional: Process the Metacell Seurat Object
由于我们将 Metacell 表达信息存储为单独的 Seurat 对象，因此我们可以对 Metacell 数据运行 Seurat 函数。 我们可以使用 GetMetacellObject 从 hdWGCNA 实验中获取元单元对象。
metacell_obj <- GetMetacellObject(seurat_obj)
Additionally, we have included a few wrapper functions to apply the Seurat workflow to the metacell object within the hdWGCNA experiment. Here we apply these wrapper functions to process the metacell object and visualize the aggregated expression profiles in two dimensions with UMAP.
seurat_obj <- NormalizeMetacells(seurat_obj)
seurat_obj <- ScaleMetacells(seurat_obj, features=VariableFeatures(seurat_obj))
seurat_obj <- RunPCAMetacells(seurat_obj, features=VariableFeatures(seurat_obj))
seurat_obj <- RunHarmonyMetacells(seurat_obj, group.by.vars='Sample')
seurat_obj <- RunUMAPMetacells(seurat_obj, reduction='harmony', dims=1:15)
这些结果储存在Seurat object@misc slot
p1 <- DimPlotMetacells(seurat_obj, group.by='cell_type') + umap_theme() + ggtitle("Cell Type")
p2 <- DimPlotMetacells(seurat_obj, group.by='Sample') + umap_theme() + ggtitle("Sample")
p1 | p2
3. 共表达网络分析
3.1 Set up the expression matrix
这一步将指定将用于网络分析的表达矩阵，使用SetDatExpr 函数，用于存储将用于下游网络分析的给定单元组的转置表达式矩阵。这一步默认情况下使用元单元表达式矩阵 (use_metacells=TRUE)。
由于我们只想分析抑制性神经元 (INH) ，因此使用group_name进行了指定。如果想包含多个细胞类型，可以用group_name = c("INH", "EX")
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
)
3.2 挑选软阈值
接下来我们将选择“软阈值”。这是 hdWGNCA 分析（以及普通 WGCNA）中极其重要的一步。 hdWGCNA 构建基因-基因相关邻接矩阵来推断基因之间的共表达关系。将相关性提高到一个幂以减少相关矩阵中存在的噪声，从而保留强连接并去除弱连接。因此，确定软功率阈值的适当值至关重要。
我们使用 TestSoftPowers 函数来为不同的软阈值进行parameter sweep。此功能通过检查不同power值的结果网络拓扑结构，帮助我们指导我们选择构建共表达网络的软阈值。共表达网络应该具有无标度拓扑，因此 TestSoftPowers 函数模拟了共表达网络在不同软功率阈值下与无标度图的相似程度。此外，我们可以使用函数 PlotSoftPowers 来可视化parameter sweep的结果。
# 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)
WGCNA 和 hdWGCNA 的一般要求是选择无标度拓扑模型拟合大于或等于 0.8 的最低软阈值，因此在这种情况下，我们选择9为软阈值。如果不设置软阈值，ConstructNetwork 将 自动选择软阈值。
Parameter sweep的输出结果存储在 hdWGCNA experiment中，可以使用 GetPowerTable 函数访问以进行进一步检查：
power_table <- GetPowerTable(seurat_obj)
head(power_table)
# Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
# 1 1 0.192964167 10.1742363 0.9508966 5887.9290 5901.7852 6501.0283
# 2 2 0.001621749 0.4801145 0.9828730 3092.6379 3092.7878 3835.2202
# 3 3 0.150812707 -3.2345892 0.9918393 1657.6198 1646.6169 2349.6967
# 4 4 0.407513393 -4.5383131 0.9898440 905.8402 890.9431 1486.7695
# 5 5 0.585949838 -4.8303680 0.9917528 504.4288 490.2683 968.1157
# 6 6 0.677668481 -4.6749295 0.9935437 286.1585 274.3233 646.6741
3.3 构建共表达网络
随后我们就可以使用ConstructNetwork函数构建共表达网络了。（The parameters for blockwiseConsensusModules can be passed directly to ConstructNetwork with the same parameter names.）
# 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')
模块组成基因
seurat_obj@misc$tutorial$wgcna_modules %>% head
# gene_name module color kME_grey kME_black kME_green kME_grey60 kME_red kME_lightcyan kME_tan kME_royalblue kME_brown kME_turquoise kME_lightyellow
# AL627309.1 AL627309.1 grey grey 0.06467554 0.01169450 0.01346971 -0.03798701 0.06389461 -0.0234879179 0.024303002 0.038220179 0.04439502 0.06965388 0.032234528
# LINC01409 LINC01409 black black 0.09562911 0.14258271 0.04532450 -0.03838167 0.07259315 0.0088322467 0.089630925 0.090604216 0.05090351 0.08046802 0.020385146
# LINC01128 LINC01128 grey grey 0.15025185 0.01456062 0.08874896 0.05263830 0.10848233 0.0233505152 0.036381392 0.003267192 0.03539629 0.13056570 0.034278458
# NOC2L NOC2L grey grey 0.16386360 -0.03203850 0.03588840 0.02706705 0.09822920 0.0223949415 0.044098042 0.013722673 0.02077281 0.15718774 0.007470014
# AGRN AGRN grey grey 0.06157292 0.05102518 0.08148903 -0.01092934 0.07061436 -0.0000888185 -0.011453627 0.028133157 0.06998927 0.04722421 0.043001474
# C1orf159 C1orf159 grey grey 0.03297087 0.02592184 0.02742151 -0.01461799 0.01558220 -0.0035206063 0.001851281 0.010550009 0.01190812 0.02965008 0.025278788
# kME_lightgreen kME_yellow kME_salmon kME_pink kME_magenta kME_purple kME_cyan kME_greenyellow kME_midnightblue kME_blue kME_darkred kME_darkgreen
# AL627309.1 0.04998851 0.03999694 0.05183194 0.057363638 0.02197360 0.07381381 0.022132715 0.0309476545 0.020243486 0.039376677 0.04946240 -0.041431912
# LINC01409 0.06062908 0.12675277 0.06172935 0.086179445 0.03545920 0.10836689 -0.010795315 0.1022885505 -0.013422943 0.008179931 0.12100397 -0.005445180
# LINC01128 0.12536877 0.04878530 0.09510907 0.037100032 0.06866021 0.07657537 0.108189685 0.0433716590 0.045969133 0.112251246 0.08318723 -0.045492473
# NOC2L 0.07886350 0.01475673 0.13338651 0.049125991 0.07227435 0.07841803 0.157399769 0.0232836501 0.025693226 0.166539003 0.03711394 -0.050635030
# AGRN 0.05589207 0.11086059 0.02339966 0.032878117 -0.03134549 0.05606373 0.036350264 0.0422577304 0.009418002 0.034792567 0.01463038 0.012060664
# C1orf159 0.01596504 0.03185948 0.02135547 0.004147534 0.00511424 0.02773708 -0.003693451 0.0004414052 0.005552397 -0.002849447 0.03548955 0.002986991
table(seurat_obj@misc$tutorial$wgcna_modules$module)
# grey black green grey60 red lightcyan tan royalblue brown turquoise lightyellow lightgreen
# 9597 154 188 87 166 88 129 67 230 281 76 86
# yellow salmon pink magenta purple cyan greenyellow midnightblue blue darkred darkgreen
# 203 123 151 150 149 123 131 112 232 60 56
3.4 Optional: inspect the topoligcal overlap matrix (TOM)
hdWGCNA represents the co-expression network as a topoligcal overlap matrix (TOM). This is a square matrix of genes by genes, where each value is the topoligcal overlap between the genes. The TOM is written to the disk when running ConstructNetwork, and we can load it into R using the GetTOM function. Advanced users may wish to inspect the TOM for custom downstream analyses.
TOM <- GetTOM(seurat_obj)
4. Module Eigengenes and Connectivity
在这一部分，我们将介绍如何计算单个细胞中的模块特征基因，以及如何计算每个基因的基于特征基因的连通性。
4.1 Compute harmonized module eigengenes
Module Eigengenes (MEs) 是一种常用的指标，用于总结整个共表达模块的基因表达谱。简而言之，模块特征基因是通过对包含每个模块的基因表达矩阵的子集执行主成分分析 (PCA) 来计算的。这些 PCA 矩阵中的第一个PC 就是 ME。
hdWGCNA 包含一个函数 ModuleEigengenes 来计算单个单元格中的模块特征基因。此外，我们允许用户对 ME 应用 Harmony 批量校正，从而产生协调的模块特征基因 (hME)。以下代码使用 group.by.vars 参数执行由原始样本协调的模块特征基因计算。
# need to run ScaleData first or else harmony throws an error:
seurat_obj <- ScaleData(seurat_obj, features=VariableFeatures(seurat_obj))
# compute all MEs in the full single-cell dataset
seurat_obj <- ModuleEigengenes(
seurat_obj,
group.by.vars="Sample"
) #时间较长
ME矩阵存储为一个矩阵，其中每一行是一个细胞，每一列是一个模块。 可以使用 GetMEs 函数从 Seurat 对象中提取此矩阵，该函数默认检索 hME。
# harmonized module eigengenes:
hMEs <- GetMEs(seurat_obj) #每个细胞对于每个模块的特征值
head(hMEs)
# lightcyan darkgreen red turquoise salmon purple black darkred lightgreen royalblue tan greenyellow
# AGTCTTTGTTGATTCG-11 -3.728266 -2.410602 3.8018569 5.540111 6.002367 4.829816 5.308432 4.066452 3.087318 5.677753 6.637234 7.003169
# AGCAGCCTCCAGATCA-11 -3.242143 -2.725894 3.1256861 5.243863 5.191766 4.053458 5.157324 3.919190 4.276658 5.821838 6.023351 7.905414
# CTTAGGATCTCATTCA-11 -3.697165 -1.910883 4.2795263 5.803626 5.108230 3.798885 5.144346 3.710556 3.451550 5.056269 6.701767 8.761117
# AACTCAGAGGAATGGA-11 -3.438431 -3.141482 0.9050427 2.689078 2.172585 2.790341 7.279000 4.412214 2.103700 4.536073 5.719890 6.730259
# ACGAGGAGTCTCCACT-11 -3.752210 -2.695012 3.8369253 1.600519 2.965712 3.075344 5.888933 2.554901 2.433528 5.539103 8.523871 7.847494
# ACTTACTTCGGAGCAA-11 -3.193025 -1.252224 3.3207978 4.391861 4.744844 2.290026 4.946544 2.965951 4.713082 4.360662 6.201866 6.620763
# midnightblue pink magenta brown lightyellow yellow grey green grey60 cyan blue
# AGTCTTTGTTGATTCG-11 6.576508 5.369716 6.578784 7.616044 4.476879 7.459703 37.92130 7.224577 4.975482 13.5867856 14.6903483
# AGCAGCCTCCAGATCA-11 5.517671 4.143108 5.390045 8.251945 4.413418 7.395514 38.43114 8.578032 6.748117 13.6829303 12.4873768
# CTTAGGATCTCATTCA-11 5.156724 5.510369 4.768236 8.243442 4.932943 8.103645 36.48388 7.867925 7.018552 12.0219104 10.7403398
# AACTCAGAGGAATGGA-11 5.603454 6.998764 6.598624 9.449461 6.600172 7.963021 29.11238 8.023318 1.401920 0.4015516 0.9184507
# ACGAGGAGTCTCCACT-11 7.024574 7.239748 7.934361 9.010840 5.841042 8.041680 32.85802 7.300763 4.790319 4.6499158 3.8687891
# ACTTACTTCGGAGCAA-11 5.045800 4.115106 5.537796 8.347864 6.416239 6.914440 30.72694 7.490030 5.289845 9.2345490 10.2069815
# module eigengenes:
MEs <- GetMEs(seurat_obj, harmonized=FALSE)
4.2 Compute module connectivity
在共表达网络分析中，我们经常希望关注“hub gene”，即在每个模块内高度连接的那些基因。 因此，我们希望确定每个基因的eigengene-based connectivity，也称为 kME。 hdWGCNA 使用 ModuleConnectivity 函数来计算完整单细胞数据集中的 kME 值，而不是元细胞数据集。 该函数本质上计算基因和模块特征基因之间的成对相关性。 kME可被用于该数据集中的所有细胞，但我们建议在用于计算 ConstructNetwork 的细胞类型或分组中计算 kME。
# compute eigengene-based connectivity (kME):
seurat_obj <- ModuleConnectivity(
seurat_obj,
group.by = 'cell_type', group_name = 'INH'
)
We can visualize the genes in each module ranked by kME using the PlotKMEs function.
# plot genes ranked by kME for each module
p <- PlotKMEs(seurat_obj, ncol=5)
p
4.3 Getting the module assignment table
hdWGCNA 允许使用 GetModules 函数轻松访问module assignment table。 该表由三列组成：gene_name: 存储基因的符号或 ID，module: 存储基因的模块分配，color: 存储每个模块的颜色映射，用于许多下游绘图步骤。 如果在此 hdWGCNA 实验中调用了 ModuleConnectivity，则此表将具有每个模块的 kME 的附加列。
# get the module assignment table:
modules <- GetModules(seurat_obj)
# show the first 6 columns:
head(modules[,1:6])
# gene_name module color kME_INH-M1 kME_INH-M2 kME_INH-M3
# AL627309.1 AL627309.1 grey grey -0.032349090 0.029917426 0.0379323320
# LINC01409 LINC01409 INH-M19 brown -0.045140924 -0.013473381 0.0102194168
# LINC01128 LINC01128 grey grey 0.050793988 0.109578251 0.1153173093
# NOC2L NOC2L grey grey 0.032490535 0.164524557 0.1699131451
# AGRN AGRN INH-M19 brown -0.008488577 0.035558532 0.0379326966
# C1orf159 C1orf159 grey grey -0.015618737 0.002235229 -0.0003824554
鉴定模块内hub基因
# get hub genes
hub_df <- GetHubGenes(seurat_obj, n_hubs = 10)
head(hub_df)
# gene_name module kME
# 1 SEPTIN2 INH-M1 0.3019739
# 2 MICU2 INH-M1 0.3140922
# 3 WDR37 INH-M1 0.3204428
# 4 UBE2Q2L INH-M1 0.3297472
# 5 HSD17B4 INH-M1 0.3399253
# 6 ASB3 INH-M1 0.3417400
This wraps up the critical analysis steps for hdWGCNA, so remember to save your output.
saveRDS(seurat_obj, file='hdWGCNA_object.rds')
- 计算每个细胞对于每个模块hub基因的表达活性(module score)
可使用seurat的AddModuleScore 函数或者Ucell算法
# compute gene scoring for the top 25 hub genes by kME for each module
# with Seurat method
seurat_obj <- ModuleExprScore(
seurat_obj,
n_genes = 25,
method='Seurat'
)
# compute gene scoring for the top 25 hub genes by kME for each module
# with UCell method
library(UCell)
seurat_obj <- ModuleExprScore(
seurat_obj,
n_genes = 25,
method='UCell'
)
5. Basic Visualization
5.1 Module Feature Plots
hdWGCNA includes the ModuleFeaturePlot function to consruct FeaturePlots for each co-expression module colored by each module’s uniquely assigned color.
# make a featureplot of hMEs for each module
plot_list <- ModuleFeaturePlot(
seurat_obj,
features='hMEs', # plot the hMEs
order=TRUE # order so the points with highest hMEs are on top
)
# stitch together with patchwork
wrap_plots(plot_list, ncol=6)
We can also plot the hub gene signature score using the same function:
# make a featureplot of hub scores for each module
plot_list <- ModuleFeaturePlot(
seurat_obj,
features='scores', # plot the hub gene scores
order='shuffle', # order so cells are shuffled
ucell = TRUE # depending on Seurat vs UCell for gene scoring
)
# stitch together with patchwork
wrap_plots(plot_list, ncol=6)
5.2 Module Correlations
hdWGCNA includes the ModuleCorrelogram function to visualize the correlation between each module based on their hMEs, MEs, or hub gene scores using the R package corrplot.
ModuleCorrelogram(seurat_obj)
5.3 Seurat plotting functions
The base Seurat plotting functions are also great for visualizing hdWGCNA outputs. Here we demonstrate plotting hMEs using DotPlot and VlnPlot. The key to using Seurat’s plotting functions to visualize the hdWGCNA data is to add it into the Seurat object’s @meta.data slot:
# get hMEs from seurat object
MEs <- GetMEs(seurat_obj, harmonized=TRUE)
mods <- colnames(MEs); mods <- mods[mods != 'grey']
# add hMEs to Seurat meta-data:
seurat_obj@meta.data <- cbind(seurat_obj@meta.data, MEs)
Now we can easily use Seurat’s DotPlot function:
# plot with Seurat's DotPlot function
p <- DotPlot(seurat_obj, features=mods, group.by = 'cell_type')
# 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')
# plot output
p
# plot with Seurat's DotPlot function
p <- DotPlot(seurat_obj, features=mods, group.by = 'cell_type')
# 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')
# plot output
p
6. 相关文献
参考：https://smorabit.github.io/hdWGCNA/articles/basic_tutorial.html
