我在前面也分享过几篇关于WGCNA的帖子，其实WGCNA本身就是对基因进行了加权共表达，合理分组的结果。
如果对于一个表达量矩阵，本身就至于2个分组，WT和处理，那么简单的差异分析，统计检验，我们其实就可以拿到上调或者下调的基因，各自去做富集分析。这种方式本身就可以对基因分组，并没有太多需要WGCNA的必要性。并且绝大多数样品的数量都是小于15个的，也不推荐使用WGCNA。有时候审稿的时候，也发现这好像都成为一个组学分析的常规的套路，比如简单的control或者处理，比如样品数目不多，也非要去做个WGCNA，实际上，个人觉得非常没有必要。
但是如果一个表达量矩阵里面的分组是比较多的，比如时间序列或者梯度变化等。如果每个都两两比较，那么分组就会非常多，造成上调或者下调的基因也很多，结果也会很多，这时候也难以从中发现比较interesting的生物学结论，这时候可以选择用WGCNA来试试。
也就是说，只要是多个分组，设计到多次差异分析，样品量多。这样的话，不同基因之间可以计算合理的相关性，可以根据基因的相关性划分为不同的模块了。
对于单细胞数据来说，单细胞表达量矩阵理论上每个细胞就是一列，所以那怕是单个10X样品也是好几千个细胞，如果样品数再增多，那么就是数万或者几十万细胞数量也是常见的。所以理论上可以使用单细胞表达量矩阵进行WGCNA分析，但是由于单细胞矩阵的稀疏性，需要对WGCNA方法进行调整，这个就是high dimensional WGCNA(hdWGCNA)，也成为scWGCNA做的东西。
==========安装==========
# 新建一个conda的环境
conda create -n hdWGCNA -c conda-forge r-base r-essentials
conda activate hdWGCNA
install.packages(c("Seurat", "WGCNA", "igraph", "devtools"))
devtools::install_github('smorabit/hdWGCNA', ref='dev')
=========测试数据==========
下载地址：
https://drive.google.com/drive/folders/1yxolklYrwFB9Snwr2Dp_W2eunBxaol4A
#加载所有的包
library(Seurat)
library(tidyverse)
library(cowplot)
library(patchwork)
library(WGCNA)
library(hdWGCNA)
# using the cowplot theme for ggplot
theme_set(theme_cowplot())
# set random seed for reproducibility
set.seed(12345)
seurat_obj <- readRDS('Zhou_2020.rds')
#发现后面报错了，没有做PCA，所以我补做了一个PCA
seurat_obj <- NormalizeData(seurat_obj, normalization.method = "LogNormalize", scale.factor = 1e4)
seurat_obj <- FindVariableFeatures(seurat_obj)
all.genes <- rownames(seurat_obj)
seurat_obj <- ScaleData(seurat_obj, features = all.genes)
seurat_obj<- RunPCA(seurat_obj, features = VariableFeatures(object = seurat_obj))
p<-DimPlot(seurat_obj,group.by='cell_type',label=TRUE)+umap_theme()+ggtitle('Zhou et al Control Cortex')+NoLegend()
p
先看下UMAP图的结果，以检查我们是否正确加载了数据，并确保我们已将细肌肉分为颗粒和细肌肉类型。
========运行WGCNA=======
1. 创建seurat的对象
hdWGCNA之前，需要创建seurat对象，通常存储在对象的@misc中。并且可以储存多个hdWGCNA结果，就像前面介绍seurat时候可存储不同的res的聚类结果一样，这里同样可以存储多个hdWGCNA结果。
SetupForWGCNA函数来创建seurat对象：
variable: 使用存储在修拉对象的VariableFeatures。
fraction：在整个数据集中或每个细胞中使用在特定部分细胞中表达的原因，由指定group.by。
custom：使用自定义列表表示中指定的原因。
例如：在下面这个例子中，我们将选择在该数据集中至少 5% 的细胞中表现出的原因，并将我们的 hdWGCNA 实验命名为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
)
从结果可以看出，seurat中多个misc，从中挑出了12639个基因进行后续的分析。
2. 创建metacells
创建万seurat对象后，下一步就是构建metacells。单细胞数据具有严重的稀疏性（data sparsity）：每个细胞中绝大部分基因的表达值为0，该比例甚至可以高达90%。但是WGCNA分析对数据的稀疏性及其敏感，因此需要构建metacells，来代替原始表达矩阵。
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可以使用参数调整元单元之间允许的重叠量。
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
reduction="harmony",
ident.group = "cell_type" # set the Idents of the metacell seurat object
)
# normalize metacell expression matrix:
seurat_obj <- NormalizeMetacells(seurat_obj)
由于我们将 Metacell 表达式信息存储为它自己的 Seurat 对象，因此我们可以在元单元格数据上运行 Seurat 函数。我们可以使用从 hdWGCNA 实验中获取 metacell 对象GetMetacellObject。
metacell_obj <- GetMetacellObject(seurat_obj)
可以看到和原始的数据对比，metacell的细胞数量大幅度减少了。
此外，我们还包含了一些包装函数，用于将 Seurat 工作流程应用于 hdWGCNA 实验中的元细胞对象。在这里，我们应用这些包装函数来处理元细胞对象，并使用 UMAP 在二维中可视化聚合的表达谱。
seurat_obj <- NormalizeMetacells(seurat_obj)
seurat_obj <- FindVariableFeatures(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
3. 共表达网络分析
我们将使用 hdWGNCA 对示例数据集中的抑制性神经元 (INH) 细胞执行共表达网络分析。
3.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
)
如果是多组值的话：group_name = c("INH", "EX")。
3.2 选择SoftPowers
这是 hdWGNCA 管道（也是普通 WGCNA）中极其重要的一步。hdWGCNA 构建了一个基因-基因相关邻接矩阵来推断基因之间的共表达关系。将相关性提高到幂以减少相关矩阵中存在的噪声量，从而保留强连接并移除弱连接。因此，确定合适的软功率阈值非常重要。
我们包含一个函数TestSoftPowers，用于针对不同的软功率阈值执行参数扫描。此功能帮助我们通过检查不同功率值的结果网络拓扑来指导我们选择构建共表达网络的软功率阈值。共表达网络应具有无标度拓扑结构，因此该TestSoftPowers函数模拟共表达网络在不同软功率阈值下与无标度图的相似程度。此外，我们还包括一个函数PlotSoftPowers来可视化参数扫描的结果。
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 的最低软功率阈值，因此在这种情况下，我们将选择我们的软功率阈值 8。如果不设置阈值，ConstructNetwork将自动选择。
参数扫描的输出表存储在 hdWGCNA 实验中，可以使用GetPowerTable函数访问以进行进一步检查：
power_table <- GetPowerTable(seurat_obj)
head(power_table)
也可以自己从表中查找。
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]
3.3 构建共表达网络
我们现在拥有构建共表达网络所需的一切。这里我们使用 hdWGCNA 函数ConstructNetwork，它在后台调用 WGCNA 函数blockwiseConsensusModules。
seurat_obj <- ConstructNetwork(
seurat_obj, soft_power=8,
setDatExpr=FALSE,
tom_name = 'INH' # name of the topoligical overlap matrix written to disk
)
hdWGCNA 还包括PlotDendrogram可视化 WGCNA 树状图的功能，这是一种常见的可视化，用于显示网络分析产生的不同共表达模块。树状图上的每片叶子代表一个基因，底部的颜色表示共表达模块分配。
PlotDendrogram(seurat_obj, main='INH hdWGCNA Dendrogram')
3.4 检查拓扑重叠矩阵 (TOM)
hdWGCNA 将共表达网络表示为拓扑重叠矩阵 (TOM)。这是一个由基因组成的方阵，其中每个值都是基因之间的拓扑重叠。有些用户可能希望返回TOM来进行下游分析，这时候就可以通过GetTOM函数来实现。
4. 模块特征和连通性
4.1 计算模块特征基因
模块特征基因 (ME) 是一种常用的指标，用于总结整个共表达模块的基因表达谱。简而言之，通过对包含每个模块的基因表达矩阵的子集执行主成分分析 (PCA) 来计算模块特征基因。每个 PCA 矩阵的第一个 PC 是 ME。
hdWGCNA 包含一个函数ModuleEigengenes来计算单个细胞中的模块特征基因。此外，我们允许用户对 ME 应用 Harmony 批量校正，从而产生协调模块特征基因 (hME)。
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)
4.2 计算模块链接
在共表达网络分析中，我们通常希望关注“中枢基因”，即在每个模块内高度连接的基因。因此，我们希望确定每个基因的基于特征基因的连通性，也称为kME。hdWGCNA 包括ModuleConnectivity计算完整单细胞数据集中的kME值，而不是元细胞数据集。这个函数本质上是计算基因和模块特征基因之间的成对相关性。
seurat_obj <- ModuleConnectivity(
seurat_obj,
group.by = 'cell_type', group_name = 'INH'
)
为方便起见，我们重新命名 hdWGCNA 模块以表明它们来自抑制性神经元组。
seurat_obj <- ResetModuleNames(
seurat_obj,
new_name = "INH-M"
)
可以使用函数可视化 kME 排序的每个模块中的基因PlotKMEs。
p <- PlotKMEs(seurat_obj, ncol=5)
p
4.3 获取模块分配表
hdWGCNA 允许使用GetModules函数轻松访问模块分配表。该表由三列组成：gene_name存储基因的符号或 ID，module存储基因的模块分配，并color存储每个模块的颜色映射，用于许多下游绘图步骤。如果ModuleConnectivity已在此 hdWGCNA 实验中调用，则此表将包含每个模块的kME的附加列。
modules <- GetModules(seurat_obj)
# get hub genes
hub_df <- GetHubGenes(seurat_obj, n_hubs = 10)
head(hub_df)
saveRDS(seurat_obj, file='hdWGCNA_object.rds')
4.4 计算hub gene的score
基因评分分析是单细胞转录组学中一种流行的方法，用于计算一组基因的整体特征的分数。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'
)
5 基本可视化
5.1 模块特征图
FeaturePlot是一种常用的 Seurat 可视化，用于直接在降维上显示感兴趣的特征。hdWGCNA 包括ModuleFeaturePlot为每个共表达模块构建 FeaturePlots 的功能，这些模块由每个模块的唯一分配颜色着色。
# make a featureplot of hMEs for each module
plot_list <- ModuleFeaturePlot(
seurat_obj,
features='hMEs', # plot the hMEs 这里可以展示四种分数hMEs,MEs,scores, average
order=TRUE # order so the points with highest hMEs are on top
)
# stitch together with patchwork
wrap_plots(plot_list, ncol=6)
5.2 模块相关性
hdWGCNA 包括使用 R 包corrplotModuleCorrelogram可视化每个模块之间基于它们的 hME、ME 或 hub 基因评分的相关性的功能。
# plot module correlagram
ModuleCorrelogram(seurat_obj)
5.3 气泡图
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)
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
5.4 小提琴图
# 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
也可以批量出图：
plot_list <- lapply(mods, function(x){
print(x)
p <- VlnPlot(
seurat_obj,
features = x,
group.by = 'cell_type',
pt.size = 0
)
p <- p + geom_boxplot(width=.25, fill='white')
p <- p + xlab('') + ylab('hME') + NoLegend()
})
wrap_plots(plot_list, ncol=6)
6. 网络可视化
6.1 单个模块网络图
library(igraph)
ModuleNetworkPlot(
seurat_obj,
mods="all",
outdir="ModuleNetworks",
plot_size=c(6,6),
label_center=FALSE,
edge.alpha=0.25,
vertex.label.cex=1,
vertex.size=6
)
输出文件夹中会有每个module的网络图。kMED的前十个基因位于图的中心，其余基因放置在外圈。
6.2 组合网络图
HubGeneNetworkPlot(
seurat_obj,
mods="all",
n_hubs=3,
n_other=6,
edge_prop=0.75,
)
这样可以可视化网络中的前3个中心基因以及6个其它基因。
6.3 UMAP 网络图
另一种方法是UMAP来可视化共表达网络中所有的基因，主要在topological overlap matrix(TOM)矩阵来计算UMAP坐标：
seurat_obj <- RunModuleUMAP(
seurat_obj,
n_hubs=10,
n_neighbors=15,
min_dist=0.1
)
umap_df <- GetModuleUMAP(seurat_obj)
ggplot(umap_df,aes(x=UMAP1,y=UMAP2))+
geom_point(color=umap_df$color,size=umap_df$kME*2)+
umap_theme()
ModuleUMAPPlot(
seurat_obj,
edge.alpha=0.25,
sample_edges=TRUE,
edge_prop=0.1,
label_hubs=2,
keep_grey_edges=FALSE
)
7. 模块与性状的相关性
seurat_obj$msex <- as.factor(seurat_obj$msex)
seurat_obj$age_death <- as.numeric(seurat_obj$age_death)
cur_traits <- c('braaksc','pmi','msex','age_death','doublet_scores','nCount_RNA','nFeature_RNA','total_counts_mt')
str(seurat_obj@meta.data[,cur_traits])
seurat_obj <- ModuleTraitCorrelation(
seurat_obj,
traits=cur_traits,
features="hMEs",
cor_method="pearson",
group.by='cell_type'
)
mt_cor <- GetModuleTraitCorrelation(seurat_obj)
names(mt_cor)
PlotModuleTraitCorrelation(
seurat_obj,
label='fdr',
label_symbol='stars',
text_size=2,
text_digits=2,
text_color='white',
high_color='#fc9272',
mid_color='#ffffbf',
low_color='#9ecae1',
plot_max=0.2,
combine=T
)
PlotModuleTraitCorrelation(
seurat_obj,
label='fdr',
label_symbol='stars',
text_size=2,
text_digits=2,
text_color='white',
high_color='#fc9272',
mid_color='#ffffbf',
low_color='#9ecae1',
plot_max=0.2,
combine=F
)
这个图实在是难看，也可以把数据导出，自己画图。
