作者，Evil Genius
偷个懒，不看文章了，总结一下关于WGCNA的分析内容。
WGCNA运用很广，在文章中一般作为补充分析存在，例如下面的结果
还有单细胞空间联合的WGCNA分析
我们就来梳理一下有关WGCNA（包括scWGCNA & hdWGNCA）在单细胞空间数据中的运用。
WGCNA（Weighted Gene Co-expression Network Analysis）是一种强大的生物信息学工具，用于分析基因表达数据，特别是识别基因间的共表达模式，寻找与特定表型（如疾病、临床特征）相关的基因模块，并提供对生物学过程的深入理解。其主要目标是通过构建加权基因共表达网络，找出在不同样本或条件下表现出相似表达模式的基因模块，进而揭示基因与生物学现象的关系。
最开始WGCNA是为了bulk-seq准备的，主要的分析内容包括：
1、基因模块的识别
2、与临床表型的关联分析
3、生物学功能分析
4、多组学数据整合
5、疾病相关基因的筛选
6、单细胞数据的应用
1. 基因模块的识别
WGCNA 最重要的应用之一是 基因模块的识别，即识别在不同样本（如不同疾病状态、不同时间点等）下表现出相似表达模式的基因集合。基因模块通常反映了某些特定的生物学过程或功能通路。
步骤：
计算基因相关性矩阵：首先计算每对基因之间的相关性，通常使用 Pearson 相关系数。
构建加权网络：使用相关性矩阵构建加权网络（每对基因之间的边权重与基因间的相关性成正比）。
识别基因模块：通过层次聚类分析识别基因模块，即将在网络中密切连接的基因聚集到一起，形成模块。
应用案例：
癌症研究：通过分析癌症样本中的基因共表达模式，识别与癌症相关的基因模块。例如，可以通过 WGCNA 找到与肿瘤侵袭性、免疫逃逸或化疗耐药性相关的基因模块。
植物基因组学：识别在不同环境条件下调控植物生长的基因模块，例如水分应激响应、光照调控等。
2. 与临床表型的关联分析
WGCNA 的另一重要应用是通过分析基因模块与 临床表型（如疾病状态、临床特征、患者生存期等）的关系，揭示基因模块与疾病的关联。这一方法常用于 疾病生物标志物 的发现。
步骤：
计算模块-表型关联：通过将每个基因模块的 模块特征值（即模块中基因的加权平均表达值）与临床数据进行相关分析，评估每个基因模块与临床表型的相关性。
模块选择与富集分析：根据模块与表型的相关性进行筛选，选出与临床特征显著相关的模块，并进行 GO 和 KEGG 富集分析，探究这些模块的生物学意义。
应用案例：
癌症生物学：例如，在乳腺癌、肺癌等癌症研究中，可以通过 WGCNA 识别与患者生存期、免疫反应、肿瘤分期等临床特征相关的基因模块。
神经退行性疾病：在阿尔茨海默病或帕金森病的研究中，WGCNA 被用于识别与认知功能、病理特征（如神经纤维缠结、β-淀粉样蛋白积累）等表型相关的基因模块。
3. 生物学功能分析
一旦识别出基因模块，可以对这些模块进行 功能富集分析，揭示它们在特定生物学过程中所起的作用。这通常通过 GO（Gene Ontology） 和 KEGG（Kyoto Encyclopedia of Genes and Genomes） 等数据库进行。
步骤：
功能富集分析：对于与表型或疾病相关的模块，使用 GO 和 KEGG 富集分析工具，推测这些模块可能涉及的生物学过程、分子功能、细胞组分以及相关的信号通路。
应用案例：
肿瘤免疫：通过对与肿瘤免疫反应相关的基因模块进行功能富集分析，揭示可能与肿瘤免疫逃逸相关的基因通路。
细胞周期调控：在研究细胞增殖和分裂的过程中，WGCNA 可以帮助识别与细胞周期相关的基因模块，并通过富集分析找到这些基因涉及的调控机制。
4. 多组学数据整合
WGCNA 还可以用于 多组学数据整合，将基因表达数据与其他组学数据（如 表观遗传学数据、蛋白质组学数据、代谢组学数据）结合，帮助揭示不同层次的生物学调控网络。
步骤：
跨组学数据整合：将不同组学的数据整合后，构建基因共表达网络，找出不同数据层次之间的共表达模块。
跨组学功能富集分析：根据不同组学的数据，分析模块的功能，并对数据进行整合和验证。
应用案例：
癌症研究：在癌症的多组学分析中，WGCNA 可以整合基因表达数据与突变、表观遗传变异等数据，发现与癌症相关的调控模块。
代谢疾病：在糖尿病、肥胖症等代谢疾病的研究中，结合基因表达和代谢数据，可以发现新的疾病标志物和调控机制。
5. 疾病相关基因的筛选
WGCNA 可以帮助筛选出与疾病或临床表型高度相关的基因模块，并进一步通过 基因筛选、功能验证 等方法，发现潜在的疾病相关基因。
步骤：
筛选疾病相关模块：根据 WGCNA 识别的基因模块，与疾病或临床特征进行相关性分析，筛选出与疾病显著相关的模块。
基因筛选与验证：从相关模块中筛选出关键基因，并通过实验验证它们在疾病中的作用。
应用案例：
基因标志物发现：通过 WGCNA 识别与肿瘤患者生存期相关的基因模块，从中筛选出关键基因并进行实验验证，作为新的癌症生物标志物。
自闭症谱系障碍：通过 WGCNA 分析大脑组织的基因表达数据，筛选出与自闭症谱系障碍相关的基因模块。
6. 单细胞数据中的 WGCNA 应用
在单细胞 RNA-seq（scRNA-seq）数据中，WGCNA 可以帮助识别 细胞类型特异性 的基因模块，并进一步探讨这些模块与细胞状态或疾病的关系。
步骤：
单细胞 WGCNA：对于每个细胞（或每个细胞群体），识别在不同细胞类型中共表达的基因模块。
细胞类型与模块关联分析：通过将模块与细胞类型、发育状态、疾病等表型进行关联，找出与细胞特征相关的基因模块。
应用案例：
肿瘤免疫细胞群体：通过在肿瘤微环境中的单细胞数据中应用 WGCNA，识别与免疫细胞反应相关的基因模块。
神经发育研究：在大脑发育过程中，WGCNA 可以识别不同发育阶段的细胞类型特异性基因模块，揭示神经元发育的分子机制。
WGCNA在单细胞数据中的运用
在单细胞数据中，由于细胞数目和基因数目庞大，WGCNA 需要处理稀疏、高维数据，并能够区分细胞间的表达异质性。scWGCNA的核心目标是构建细胞间基因的加权共表达网络，并识别基因模块。
WGCNA的应用步骤
1、数据准备与预处理
- 选择感兴趣的基因和细胞，进行数据标准化和预处理。通常使用 log2 转换 和 归一化 方法来去除文库效应和技术偏差。
- 过滤低表达基因和低质量细胞。
- 对数据进行降维（如PCA）以提高计算效率。
2、构建加权基因共表达网络
- 计算每对基因之间的相关性，常使用 Pearson 相关系数。
- 使用软阈值（soft-thresholding）方法来计算加权相关性矩阵，以保证网络的平滑性。
3、模块识别与功能注释
- 对加权基因共表达网络进行 层次聚类，从而识别出具有相似表达模式的基因模块。
- 使用 GO（Gene Ontology）和 KEGG（Kyoto Encyclopedia of Genes and Genomes）等数据库进行功能富集分析，揭示模块可能的生物学意义。
4、与表型、细胞类型的关联分析
- 根据模块特征值与细胞类型或表型数据（如疾病状态、发育阶段等）进行相关性分析。
WGCNA示例代码
以下展示如何在单细胞数据中应用 WGCNA。
library(WGCNA)
library(ggplot2)
library(dplyr)
1. 数据预处理
单细胞RNA-seq的表达矩阵（expr_matrix），其中行是基因，列是细胞。需要进行数据预处理和标准化。
# 假设 expr_matrix 是一个基因-细胞的矩阵，行是基因，列是细胞
# 对数据进行log2转化（假设原始数据是count数据）
expr_matrix <- log2(expr_matrix + 1)
# 选择表达量较高的基因，去除低表达基因（比如选择平均表达值大于1的基因）
gene_filter <- rowMeans(expr_matrix) > 1
expr_matrix <- expr_matrix[gene_filter, ]
# 对数据进行标准化（例如，按细胞归一化）
expr_matrix <- t(scale(t(expr_matrix))) # 每个基因的表达值按细胞归一化
如果单细胞数据直接处理成了rds文件，则直接读取即可
data = readRDS('scRNA.rds')
expr_matrix = data[["RNA"]]@scale.data
2. 构建加权基因共表达网络
接下来，使用 WGCNA 包构建基因共表达网络，并使用软阈值选择合适的加权参数。
# 计算基因相关性矩阵（Pearson 相关系数）
correlation_matrix <- cor(expr_matrix)
# 选择合适的软阈值（一般从1开始试不同的值）
soft_power <- 6 # 假设选择了6作为软阈值
adjacency_matrix <- adjacency(correlation_matrix, power = soft_power)
# 将邻接矩阵转化为拓扑重叠矩阵（TOM）
TOM_matrix <- TOMsimilarity(adjacency_matrix)
# 通过层次聚类，识别基因模块
gene_tree <- hclust(as.dist(1 - TOM_matrix), method = "average")
plot(gene_tree, main = "Gene Dendrogram", xlab = "", sub = "")
# 基于树状图剪枝，识别基因模块
dynamic_modules <- cutreeDynamic(dendro = gene_tree, distM = 1 - TOM_matrix, deepSplit = 2, pamRespectsDendro = FALSE)
table(dynamic_modules)
3. 模块特征值的计算
计算每个基因模块的模块特征值（即模块中基因的加权平均表达值）。
# 计算模块特征值
module_trait_relationship <- moduleTraitCor(cor(expr_matrix, dynamic_modules))
4. 与临床表型或细胞类型的关联分析(单细胞主要是和细胞类型相关联，尤其是细胞亚型)
假设你有一个细胞类型或临床表型的数据，接下来可以将基因模块与这些表型进行关联分析。
# 假设表型数据存储在 phenotype_data 中
# phenotype_data 可能包括细胞类型、疾病状态等信息
# 计算模块特征值与表型的相关性
module_trait_cor <- cor(module_trait_relationship, phenotype_data)
# 可视化关联结果
heatmap(module_trait_cor, main = "Module-Phenotype Relationships")
5. 功能富集分析
对识别的基因模块进行功能富集分析，了解每个模块的生物学意义。
# 假设你使用GO富集分析
library(clusterProfiler)
# 获取某个模块的基因ID
module_genes <- names(dynamic_modules)[dynamic_modules == 1] # 模块1的基因
# 进行GO富集分析
go_results <- enrichGO(gene = module_genes, OrgDb = org.Hs.eg.db, keyType = "SYMBOL", ont = "BP", pvalueCutoff = 0.05)
head(go_results)
6. 结果可视化
通过可视化工具将结果呈现，帮助分析和解释。
# 可视化模块的基因共表达网络
plotDendroAndColors(gene_tree, dynamic_modules, main = "Gene Dendrogram with Module Colors")
# 进一步可视化与临床表型相关的模块
heatmap(module_trait_cor, main = "Module-Phenotype Relationships", scale = "row")
应用案例
癌症研究：通过分析肿瘤细胞与正常细胞的基因共表达网络，识别肿瘤特有的基因模块，并与临床特征（如肿瘤分期、免疫反应等）进行关联分析。
神经发育：在大脑发育的单细胞数据中，识别不同发育阶段的细胞类型特异性基因模块，探索神经发育的分子机制。
免疫学研究：分析免疫细胞中的基因共表达模式，识别与免疫反应或免疫逃逸相关的基因模块。
scWGCNA在单细胞数据中的运用
scWGCNA 是一个专门为单细胞数据设计的加权基因共表达网络分析方法。它通过考虑单细胞数据的特点，如细胞异质性、基因表达的高变性等，来建立基因共表达网络并识别具有共表达关系的基因模块。在单细胞RNA-seq分析中，scWGCNA可以帮助揭示基因模块的空间或细胞类型特异性，从而深入了解基因功能和调控机制。
# 安装devtools（如果没有安装的话）
#install.packages("devtools")
# 使用devtools从GitHub安装scWGCNA
#devtools::install_github("YosefLab/scWGCNA")
# 加载相关包
library(scWGCNA)
library(Seurat)
library(clusterProfiler)
library(org.Hs.eg.db) # 如果是人类基因数据
# 加载Seurat对象
seurat_obj <- readRDS("path_to_your_seurat_object.rds")
# 提取标准化后的表达矩阵
expr_matrix <- seurat_obj@assays$RNA@data
# 获取基因和细胞类型的元数据
cell_types <- seurat_obj@meta.data$cell_type # 假设细胞类型存储在meta.data中
####构建加权基因共表达网络
# 选择高变基因（默认选择前2000个高变基因）
high_var_genes <- VariableFeatures(seurat_obj)
# 提取这些高变基因的表达矩阵
expr_matrix_high_var <- expr_matrix[high_var_genes, ]
# 使用scWGCNA的函数构建加权基因共表达网络
network <- scWGCNA::buildNetwork(expr_matrix_high_var)
###识别基因模块
# 识别基因模块
modules <- scWGCNA::findModules(network)
# 查看识别的模块信息
modules$moduleAssignments # 每个基因所属的模块
###计算模块特征基因值 (Module Eigengene, ME)
# 计算模块特征基因值
module_eigengenes <- scWGCNA::calculateEigengenes(expr_matrix_high_var, modules$moduleAssignments)
# 查看模块特征基因值
head(module_eigengenes)
###进行表型分析
# 假设表型数据存储在Seurat对象的meta.data中
# 例如：细胞类型信息存储在 'cell_type' 中
cell_types <- seurat_obj@meta.data$cell_type
# 计算模块特征基因值与细胞类型的相关性
module_trait_relationship <- cor(module_eigengenes, as.numeric(cell_types), method = "pearson")
# 可视化模块与表型（例如，细胞类型）的关系
heatmap(module_trait_relationship, main = "Module vs Cell Type")
####富集分析
# 假设我们分析模块1
module_1_genes <- colnames(expr_matrix)[modules$moduleAssignments == 1]
# 使用clusterProfiler进行GO富集分析
go_results <- enrichGO(
gene = module_1_genes,
OrgDb = org.Hs.eg.db, # 选择适当的物种数据库
keyType = "SYMBOL", # 基因符号（如果是ENSEMBL ID等其他ID类型，替换）
ont = "BP", # 生物过程（"BP"）/细胞组件（"CC"）/分子功能（"MF"）
pAdjustMethod = "BH", # 多重检验校正
qvalueCutoff = 0.05 # 设置q值阈值
)
# 查看GO富集结果
summary(go_results)
# 可视化GO富集结果
barplot(go_results, showCategory = 10)
###可视化模块表达模式
# 绘制模块特征基因值在不同细胞类型中的箱型图
library(ggplot2)
# 假设我们选择模块1
module_1_ME <- module_eigengenes[, 1]
# 创建数据框
module_df <- data.frame(
cell_type = cell_types,
module_1_ME = module_1_ME
)
# 绘制箱型图
ggplot(module_df, aes(x = cell_type, y = module_1_ME)) +
geom_boxplot() +
labs(title = "Module 1 ME across Cell Types")
hdWGCNA在单细胞空间数据中的运用
hdWGCNA 是一种扩展的加权基因共表达网络分析（WGCNA）方法，专门用于处理单细胞空间数据（例如空间转录组学数据）。与传统的WGCNA不同，hdWGCNA 可以结合空间信息，考虑基因在空间中的共表达模式。这使得hdWGCNA特别适合于空间转录组学数据的分析，在这些数据中，基因表达不仅受到细胞类型和状态的影响，还受到组织或区域的空间位置影响。
hdWGCNA在单细胞空间数据中的运用
在空间转录组学中，细胞的空间位置和基因表达具有重要关系。通过 hdWGCNA，可以：
- 识别空间基因模块：分析空间相关性，并将基因分为多个模块。
- 探讨空间与表型的关系：通过空间基因模块与细胞类型、疾病状态等表型数据之间的相关性分* 析，揭示空间基因模块的功能和生物学意义。
- 空间信息集成：不仅考虑基因共表达关系，还将空间信息作为网络的一部分进行分析。
分析框架
1、数据准备与加载
- 准备单细胞空间转录组数据，数据通常包含基因表达矩阵、细胞空间坐标等信息。
2、构建加权基因共表达网络- 根据基因之间的相关性构建网络，考虑基因在空间中的表达特征。
- 需要选择合适的软阈值（soft-thresholding）来建立加权网络。
3、模块识别与空间特征分析- 对加权网络进行聚类，识别基因模块。
- 将空间位置与基因模块关联，探讨基因模块在空间中的分布。
4、模块与表型关联分析- 根据基因模块特征值与表型数据进行相关性分析，探索空间模块与细胞类型或疾病等表型之间的关系。
# 安装 hdWGCNA 包
devtools::install_github("hdWGCNA/hdWGCNA")
# 加载 hdWGCNA 包
library(hdWGCNA)
# 加载 Seurat 包并读取对象
library(Seurat)
seurat_obj <- readRDS("path_to_your_seurat_object.rds")
# 获取基因表达矩阵（已经归一化）
expr_matrix <- seurat_obj@assays$RNA@data
# 获取细胞的空间坐标（假设空间坐标存储在 metadata 中）
cell_coords <- seurat_obj@meta.data[, c("spatial_dim1", "spatial_dim2")] # 适当替换列名
# 构建加权基因共表达网络
network <- hdWGCNA::createNetwork(
exprData = expr_matrix, # 基因表达矩阵
spatialData = cell_coords, # 空间坐标
power = 6 # 选择软阈值
)
# 查看网络的基本信息
summary(network)
# 识别基因模块
modules <- hdWGCNA::identifyModules(
network = network # 使用前面创建的网络
)
# 查看模块信息
print(modules)
###空间信息分析
# 提取模块特征基因值（模块的代表基因）
module_eigengenes <- hdWGCNA::getModuleEigengenes(
network = network, # 使用前面创建的网络
modules = modules # 使用识别的基因模块
)
# 可视化模块特征基因值在空间坐标上的分布
library(ggplot2)
space_df <- data.frame(
x = cell_coords[, 1],
y = cell_coords[, 2],
module_1 = module_eigengenes[, 1], # 假设你有多个模块
module_2 = module_eigengenes[, 2]
)
# 绘制模块1在空间中的分布
ggplot(space_df, aes(x = x, y = y, color = module_1)) +
geom_point() +
labs(title = "Spatial Distribution of Module 1")
####表性分析
# 假设你有一个表型数据（例如细胞类型），并将其存储在 Seurat 对象的 metadata 中
phenotype_data <- seurat_obj@meta.data$cell_type # 替换为实际的表型数据列
# 计算模块特征基因值与表型的相关性
module_trait_relationship <- cor(module_eigengenes, phenotype_data)
# 可视化模块与表型的关系
heatmap(module_trait_relationship, main = "Module-Phenotype Relationships")
####模块富集分析
# 安装和加载需要的R包
install.packages("clusterProfiler")
install.packages("org.Hs.eg.db") # 以人类基因为例
install.packages("enrichplot")
install.packages("ggplot2")
library(clusterProfiler)
library(org.Hs.eg.db)
library(enrichplot)
library(ggplot2)
# 提取模块特征基因值（代表基因）
module_eigengenes <- hdWGCNA::getModuleEigengenes(
network = network, # 使用前面创建的网络
modules = modules # 使用前面识别的模块
)
# 提取每个模块的基因集合（假设模块名为 modules）
module_genes <- lapply(modules, function(mod) {
gene_names <- colnames(expr_matrix)[which(gene_modules == mod)] # 找到该模块的基因
return(gene_names)
})
# 查看某个模块的基因
head(module_genes[[1]])
# 选择某个模块（假设我们分析模块1）
module1_genes <- module_genes[[1]]
# GO富集分析
go_results <- enrichGO(
gene = module1_genes,
OrgDb = org.Hs.eg.db, # 选择人类基因数据库
keyType = "SYMBOL", # 基因符号（如果是其他类型的ID，如ENSEMBL，可以选择对应类型）
ont = "BP", # 生物过程（"BP"）/细胞组件（"CC"）/分子功能（"MF"）
pAdjustMethod = "BH", # 多重检验校正方法
qvalueCutoff = 0.05 # 设置q值阈值
)
# 查看GO富集结果
summary(go_results)
# KEGG富集分析
kegg_results <- enrichKEGG(
gene = module1_genes,
organism = "hsa", # hsa 代表人类基因组，若是其他物种，替换为相应的代码，如 "mmu"（小鼠）
keyType = "kegg", # 使用KEGG ID（如ENTREZID）作为输入基因ID
pAdjustMethod = "BH", # 多重检验校正方法
qvalueCutoff = 0.05 # 设置q值阈值
)
# 查看KEGG富集结果
summary(kegg_results)
###可视化
# GO富集分析结果的条形图
barplot(go_results, showCategory = 10) # 显示前10个显著富集的GO项
# KEGG富集分析结果的点图
dotplot(kegg_results, showCategory = 10) # 显示前10个显著富集的KEGG通路
# 细胞过程与模块基因的网络图
cnetplot(go_results, showCategory = 5)
####综合富集分析的结果
# 对多个模块进行富集分析
go_results_all_modules <- lapply(module_genes, function(genes) {
enrichGO(gene = genes, OrgDb = org.Hs.eg.db, keyType = "SYMBOL", ont = "BP", pAdjustMethod = "BH", qvalueCutoff = 0.05)
})
# 可视化所有模块的GO富集结果
go_results_all_modules_vis <- lapply(go_results_all_modules, barplot, showCategory = 5)
# 可视化多个模块的KEGG富集结果
kegg_results_all_modules <- lapply(module_genes, function(genes) {
enrichKEGG(gene = genes, organism = "hsa", keyType = "kegg", pAdjustMethod = "BH", qvalueCutoff = 0.05)
})
