单细胞数据-hdWGCNA分析
学习整理from官网
单细胞数据中的 hdWGCNA
本教程涵盖了使用 hdWGCNA 对单细胞数据执行共表达网络分析的基础知识。在这里，我们使用来自本出版物的人类皮层样本的处理单核 RNA-seq (snRNA-seq) 数据集来演示 hdWGCNA 。该数据集已经使用标准的单细胞转录组学分析管道（例如Seurat或Scanpy ）进行了全面处理。如果您想使用自己的数据集学习本教程，您首先需要满足以下先决条件：
Seurat 格式的单细胞或单核转录组学数据集。
标准化基因表达矩阵NormalizeData。
识别高度可变的基因VariableFeatures。
缩放归一化的表达式数据ScaleData
如果需要，执行降维RunPCA和批量校正RunHarmony。
用于可视化的非线性降维RunUMAP。
将细胞分组为簇（FindNeighbors和FindClusters）。
可以在Seurat Guided Clustering Tutorial中找到运行先决条件数据处理步骤的示例。
此外，还有很多 WGCNA 特定的术语和首字母缩略词，都在该表中进行了说明。
1. 下载教程数据
出于本教程的目的，我们提供了周等人2020年研究中控人脑的经过处理的修拉对象。
wget https://swaruplab.bio.uci.edu/public_data/Zhou_2020.rds
下载不能使用？
请尝试从这里链接下载文件。
2. 添加数据集和所需的库
首先，我们将加载单细胞数据集和本教程所需的 R 库。
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)
# load the Zhou et al snRNA-seq dataset
seurat_obj <- readRDS('Zhou_2020.rds')
在这里，我们将绘制按细肌肉类型着着颜色的UMAP，以检查我们是否正确加载了数据，并确保我们已将细肌肉分为颗粒和细肌肉类型。
p <- DimPlot(seurat_obj, group.by='cell_type', label=TRUE) +
umap_theme() + ggtitle('Zhou et al Control Cortex') + NoLegend()
p
image.png
3. 为 WGCNA 设置 Seurat 对象
hdwgcna之前，seurat seurat对象hdwgcna计算计算计算的大部分信息信息都在存储存储存储存储存储存储存储对象对象对象对象对象对象的@misc槽槽槽中。值得注意的是，由于我们将hdWGCNA看为下游数据分析步骤，因此我们不支持在运行后对Seurat对象进行子集化SetupForWGCNA。
SetupForWGCNAseurat seurat对象gene_select对象
- variable: 使用存储在修拉对象的VariableFeatures。
- fraction：在整个数据集中或每个细胞中使用在特定部分细胞中表达的原因，由指定group.by。
- custom：使用自定义列表表示中指定的原因。
在这个例子中，我们将选择在该数据集中至少 5% 的细胞中表现出的原因，并将我们的 hdWGCNA 实验命名为“教程”。
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
)
4. 构建元细胞
设置好 Seurat 对象后，在 hdWGCNA 中运行 hdWGCNA 管道的第一步是从单细胞数据集中构建元细胞。简而言之，元细胞是源自相同生物来源样本的小群相似细胞的聚集体。k-最近邻 (KNN) 算法用于识别要聚集的相似细胞组，然后计算这些细胞的平均或求和表达，从而产生元细胞基因表达矩阵。与原始表达矩阵相比，元细胞表达矩阵的稀疏性大大降低，因此更适合使用。我们最初的动机是使用元单元代替原始的单个单元，因为 WGCNA 等相关网络方法对数据稀疏性很敏感。此外，Cicero，在构建共同可访问性网络之前采用类似的元单元聚合方法。
hdWGCNA 包含一个函数MetacellsByGroups，用于在给定单细胞数据集的情况下构建元细胞表达矩阵。此函数为存储在 hdWGCNA 实验内部的 metacell 数据集构造了一个新的 Seurat 对象。该group.by参数决定了将在哪些组中构建元细胞。我们只想从来自相同来源生物样本的细胞构建元细胞，因此通过参数将该信息传递给 hdWGCNA 至关重要group.by。此外，我们通常分别为每种细胞类型构建元细胞。因此，在此示例中，我们按Sample和进行分组cell_type以获得所需的结果。
要聚合的单元格数量k应根据输入数据集的大小进行调整，通常较小的数量k可用于小型数据集。我们通常使用k20 到 75 之间的值。本教程使用的数据集有 40,039 个细胞，每个生物样本的细胞数从 890 到 8,188 不等，这里我们使用k=25. max_shared可以使用参数调整元单元之间允许的重叠量。
注意：我们发现元细胞聚合方法对于代表性极低的细胞类型不会产生良好的结果。例如，在该数据集中，脑血管细胞（周细胞和内皮细胞）的代表最少，我们已将它们排除在该分析之外。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)
可选：处理 Metacell Seurat 对象
由于我们将 Metacell 表达式信息存储为它自己的 Seurat 对象，因此我们可以在元单元格数据上运行 Seurat 函数。我们可以使用从 hdWGCNA 实验中获取 metacell 对象GetMetacellObject。
metacell_obj <- GetMetacellObject(seurat_obj)
此外，我们还包含了一些包装函数，用于将 Seurat 工作流程应用于 hdWGCNA 实验中的元细胞对象。在这里，我们应用这些包装函数来处理元细胞对象，并使用 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)
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
image.png
5. 共表达网络分析
在本节中，我们将讨论如何使用 hdWGNCA 对示例数据集中的抑制性神经元 (INH) 细胞执行共表达网络分析。
5.1 建立表达矩阵
此处我们指定将用于网络分析的表达式矩阵。因为我们只想包括抑制神经元，所以我们必须在构建网络之前对我们的表达数据进行子集化。hdWGCNA 包括SetDatExpr为给定的一组细胞存储转置表达矩阵的功能，这些细胞将用于下游网络分析。默认情况下使用元细胞表达矩阵 ( use_metacells=TRUE)，但 hdWGCNA 确实允许在需要时使用单细胞表达矩阵。此功能允许用户指定从哪个插槽获取表达矩阵，例如，如果用户想要应用SCTransform规范化而不是NormalizeData.
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
)
选择多个组
假设您要同时对多种细胞类型或簇进行共表达网络分析。SetDatExpr可以使用略有不同的设置运行，通过将字符向量传递给group_name参数来实现所需的结果。
seurat_obj <- SetDatExpr(
seurat_obj,
group_name = c("INH", "EX"),
group.by='cell_type'
)
5.2 选择"软实力门槛"[SoftPowers]
选择软实力门槛
接下来我们将选择“软实力门槛”。这是 hdWGNCA 管道（以及普通 WGCNA）中极其重要的一步。hdWGCNA 构建了一个基因-基因相关邻接矩阵来推断基因之间的共表达关系。将相关性提高到幂以减少相关矩阵中存在的噪声量，从而保留强连接并移除弱连接。因此，确定合适的软功率阈值非常重要。
我们包含一个函数TestSoftPowers，用于针对不同的软功率阈值执行参数扫描。此功能帮助我们通过检查不同功率值的结果网络拓扑来指导我们选择构建共表达网络的软功率阈值。共表达网络应具有无标度拓扑结构，因此该TestSoftPowers函数模拟共表达网络在不同软功率阈值下与无标度图的相似程度。此外，我们还包括一个函数PlotSoftPowers来可视化参数扫描的结果。
以下代码执行参数扫描并输出摘要图。
# 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)
image.png
WGCNA 和 hdWGCNA 的一般指导是选择无标度拓扑模型拟合大于或等于 0.8 的最低软功率阈值，因此在这种情况下，我们将选择我们的软功率阈值 9。稍后，ConstructNetwork将自动如果用户没有提供，请选择软功率阈值。
参数扫描的输出表存储在 hdWGCNA 实验中，可以使用GetPowerTable函数访问以进行进一步检查：
power_table <- GetPowerTable(seurat_obj)
head(power_table)
输出
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
a=power_table$SFT.R.sq
i = 1
for (b in a){
print(b)
if (b<0.8) {
i=i+1
print(i)
}else if(b>0.8){
break
}
}
select_soft_power = power_table$Power[i]
5.3 构建共表达网络
我们现在拥有构建共表达网络所需的一切。这里我们使用 hdWGCNA 函数ConstructNetwork，它在后台调用 WGCNA 函数blockwiseConsensusModules。如果您是高级用户，此函数有很多参数可供使用，但我们选择了适用于许多单细胞数据集的默认参数。for 的参数可以使用相同的参数名称blockwiseConsensusModules直接传递给。ConstructNetwork
以下代码使用上面选择的软功率阈值构建共表达网络：
# 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
)
hdWGCNA 还包括PlotDendrogram可视化 WGCNA 树状图的功能，这是一种常见的可视化，用于显示网络分析产生的不同共表达模块。树状图上的每片叶子代表一个基因，底部的颜色表示共表达模块分配。
重要的是，“灰色”模块由未分组到任何共表达模块中的基因组成。所有下游分析和解释都应忽略灰色模块。
PlotDendrogram(seurat_obj, main='INH hdWGCNA Dendrogram')
image.png
5.4 可选：检查拓扑重叠矩阵 (TOM)
hdWGCNA 将共表达网络表示为拓扑重叠矩阵 (TOM)。这是一个由基因组成的方阵，其中每个值都是基因之间的拓扑重叠。TOM在运行时写入磁盘，我们可以使用该函数ConstructNetwork将其加载到R中。GetTOM高级用户可能希望检查 TOM 以进行自定义下游分析。
TOM <- GetTOM(seurat_obj)
6 模块本征和连通性
在本节中，我们将介绍如何计算单个细胞中的模块特征基因，以及如何计算每个基因的基于特征基因的连通性。
6.1 计算协调模块特征基因
模块特征基因 (ME) 是一种常用的指标，用于总结整个共表达模块的基因表达谱。简而言之，通过对包含每个模块的基因表达矩阵的子集执行主成分分析 (PCA) 来计算模块特征基因。每个 PCA 矩阵的第一个 PC 是 ME。
降维技术是单细胞基因组学中非常热门的话题。众所周知，技术人为因素会使单细胞数据集的分析变得混乱，多年来，已经出现了许多旨在减少这些人为因素影响的方法。因此，按理说，ME 也会受到这些技术工件的影响，而 hdWGCNA 试图减轻这些影响。
hdWGCNA 包含一个函数ModuleEigengenes来计算单个细胞中的模块特征基因。此外，我们允许用户对 ME 应用 Harmony 批量校正，从而产生协调模块特征基因 (hME)。以下代码使用group.by.vars参数执行由 Sample of origin 协调的模块 eigengene 计算。
# 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可以使用默认检索hME的函数从 Seurat 对象中提取该矩阵。
# harmonized module eigengenes:
hMEs <- GetMEs(seurat_obj)
# module eigengenes:
MEs <- GetMEs(seurat_obj, harmonized=FALSE)
6.2 计算模块连接
计算模块连接
在共表达网络分析中，我们通常希望关注“中枢基因”，即在每个模块内高度连接的基因。因此，我们希望确定每个基因的基于特征基因的连通性，也称为kME。hdWGCNA 包括ModuleConnectivity计算完整单细胞数据集中的kME值，而不是元细胞数据集。这个函数本质上是计算基因和模块特征基因之间的成对相关性。可以为数据集中的所有细胞计算 kME，但我们建议在以前用于运行的细胞类型或组中计算 kME ConstructNetwork。
# compute eigengene-based connectivity (kME):
seurat_obj <- ModuleConnectivity(
seurat_obj,
group.by = 'cell_type', group_name = 'INH'
)
为方便起见，我们重新命名 hdWGCNA 模块以表明它们来自抑制性神经元组。有关重命名模块的更多信息，请参阅模块自定义教程。
# rename the modules
seurat_obj <- ResetModuleNames(
seurat_obj,
new_name = "INH-M"
)
我们可以使用函数可视化 kME 排序的每个模块中的基因PlotKMEs。
# plot genes ranked by kME for each module
p <- PlotKMEs(seurat_obj, ncol=5)
p
image.png
6.3 get the module assignment table 获取模块分配表
hdWGCNA 允许使用GetModules函数轻松访问模块分配表。该表由三列组成：gene_name存储基因的符号或 ID，module存储基因的模块分配，并color存储每个模块的颜色映射，用于许多下游绘图步骤。如果ModuleConnectivity已在此 hdWGCNA 实验中调用，则此表将包含每个模块的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
可以使用该GetHubGenes函数提取按 kME 排序的前 N 个中心基因表。
# 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
hdWGCNA 的关键分析步骤到此结束，因此请记住保存您的输出。
saveRDS(seurat_obj, file='hdWGCNA_object.rds')
6.4 计算中心基因签名分数
基因评分分析是单细胞转录组学中一种流行的方法，用于计算一组基因的整体特征的分数。Seurat 使用该函数实现了他们自己的基因评分技术AddModuleScore，但也有其他方法，例如UCell。hdWGCNA 包括ModuleExprScore使用 Seurat 或 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'
)
7.基本可视化
在这里，我们展示了 hdWGCNA 的一些基本可视化功能，并演示了如何使用 Seurat 的一些内置绘图工具来可视化我们的 hdWGCNA 结果。请注意，我们有一个单独的 hdWGCNA 网络可视化教程。
7.1 模块特征图
FeaturePlot是一种常用的 Seurat 可视化，用于直接在降维上显示感兴趣的特征。hdWGCNA 包括ModuleFeaturePlot为每个共表达模块构建 FeaturePlots 的功能，这些模块由每个模块的唯一分配颜色着色。
# 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)
image.png
7.2 模块关联
hdWGCNA 包括使用 R 包corrplotModuleCorrelogram可视化每个模块之间基于它们的 hME、ME 或 hub 基因评分的相关性的功能。
# plot module correlagram
ModuleCorrelogram(seurat_obj)
image.png
7.3 修改拉绘图函数
seurat绘图函数绘图函数非常可可视化视化视化视化视化输出输出在DotPlot这里这里，我们我们演示使用使用使用和和和绘制绘制绘制绘制绘制。。。。。。VlnPlot使用使用@meta.data
#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)
现在我们可以很容易地使用修拉的DotPlot函数：
#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
image.png
这是我们使用修拉VlnPlot函数的另一个例子：
# Plot INH-M4 hME using Seurat VlnPlot function
p <- VlnPlot(
seurat_obj,
features = 'INH-M12',
group.by = 'cell_type',
pt.size = 0 # don't show actual data points
)
# add box-and-whisker plots on top:
p <- p + geom_boxplot(width=.25, fill='white')
# change axis labels and remove legend:
p <- p + xlab('') + ylab('hME') + NoLegend()
# plot output
p
image.png
下一步
在本教程中，我们介绍了在单细胞转录组学数据中执行共表达网络分析的核心功能。我们鼓励您浏览我们的其他教程以对这些 hdWGCNA 结果进行下游分析。
