基因共表达网络分析（WGCNA）¶

library(WGCNA)
options(stringsAsFactors = FALSE)
enableWGCNAThreads(nThreads = 8)

# 读取表达矩阵（行=样本，列=基因）
# 如果原始矩阵是基因×样本，需要转置
expr_data <- read.csv("expression_matrix.csv", row.names = 1)
datExpr <- as.data.frame(t(expr_data))  # 转置：行=样本，列=基因

# 读取性状/临床信息
traitData <- read.csv("clinical_traits.csv", row.names = 1)

# === 基因过滤 ===
# 方法1：去除低表达基因（至少在一定比例样本中有表达）
# 方法2：保留变异最大的前N个基因
library(matrixStats)
gene_var <- colVars(as.matrix(datExpr))
names(gene_var) <- colnames(datExpr)
top_genes <- names(sort(gene_var, decreasing = TRUE))[1:5000]  # top 5000变异基因
datExpr <- datExpr[, top_genes]

# === 样本质控 ===
# 检查缺失值
gsg <- goodSamplesGenes(datExpr, verbose = 3)
if (!gsg$allOK) {
  datExpr <- datExpr[gsg$goodSamples, gsg$goodGenes]
}

# 样本层次聚类检测离群
sampleTree <- hclust(dist(datExpr), method = "average")
plot(sampleTree, main = "Sample Dendrogram", sub = "", xlab = "", cex.lab = 1.5)

# 去除离群样本（手动设定切割高度）
abline(h = 15, col = "red")  # 根据图形调整高度
clust <- cutreeStatic(sampleTree, cutHeight = 15, minSize = 10)
keepSamples <- (clust == 1)
datExpr <- datExpr[keepSamples, ]
traitData <- traitData[keepSamples, ]

cat("Samples:", nrow(datExpr), " Genes:", ncol(datExpr), "\n")

# 选择软阈值
powers <- c(c(1:10), seq(12, 20, by = 2))

sft <- pickSoftThreshold(datExpr,
                         powerVector = powers,
                         networkType = "signed",  # 推荐signed网络
                         verbose = 5)

# 可视化
par(mfrow = c(1, 2))

# Scale-free topology fit (R²)
plot(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2],
     xlab = "Soft Threshold (power)",
     ylab = "Scale Free Topology Model Fit (R²)",
     type = "n",
     main = "Scale Independence")
text(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2],
     labels = powers, cex = 0.9, col = "red")
abline(h = 0.85, col = "red")

# Mean connectivity
plot(sft$fitIndices[, 1], sft$fitIndices[, 5],
     xlab = "Soft Threshold (power)",
     ylab = "Mean Connectivity",
     type = "n",
     main = "Mean Connectivity")
text(sft$fitIndices[, 1], sft$fitIndices[, 5],
     labels = powers, cex = 0.9, col = "red")

# 选择power（通常选R²首次>0.85的最小power）
softPower <- sft$powerEstimate
cat("Selected soft threshold:", softPower, "\n")
# 如果自动选择不理想，手动设置（如6或8）
if (is.na(softPower)) softPower <- 8

# === 一步法（推荐初学者）===
net <- blockwiseModules(datExpr,
  power = softPower,
  networkType = "signed",
  TOMType = "signed",
  minModuleSize = 30,        # 最小模块基因数
  reassignThreshold = 0,
  mergeCutHeight = 0.25,     # 模块合并阈值（越小合并越少）
  numericLabels = TRUE,
  pamRespectsDendro = FALSE,
  saveTOMs = TRUE,
  saveTOMFileBase = "TOM",
  verbose = 3,
  maxBlockSize = 20000       # 内存限制相关
)

# 模块颜色
moduleColors <- labels2colors(net$colors)
table(moduleColors)

# 可视化模块聚类树
plotDendroAndColors(net$dendrograms[[1]],
                   moduleColors[net$blockGenes[[1]]],
                   "Module Colors",
                   dendroLabels = FALSE,
                   hang = 0.03,
                   addGuide = TRUE,
                   guideHang = 0.05)

# === 分步法（更灵活）===
# 1. 计算邻接矩阵
adjacency <- adjacency(datExpr, power = softPower, type = "signed")

# 2. 计算TOM
TOM <- TOMsimilarity(adjacency, TOMType = "signed")
dissTOM <- 1 - TOM

# 3. 层次聚类
geneTree <- hclust(as.dist(dissTOM), method = "average")

# 4. Dynamic Tree Cut
dynamicMods <- cutreeDynamic(dendro = geneTree,
                              distM = dissTOM,
                              deepSplit = 2,
                              pamRespectsDendro = FALSE,
                              minClusterSize = 30)

# 5. 合并相似模块
MEList <- moduleEigengenes(datExpr, colors = labels2colors(dynamicMods))
MEs <- MEList$eigengenes
MEDiss <- 1 - cor(MEs)
METree <- hclust(as.dist(MEDiss), method = "average")
plot(METree, main = "Clustering of module eigengenes")
abline(h = 0.25, col = "red")

merge <- mergeCloseModules(datExpr, labels2colors(dynamicMods), cutHeight = 0.25)
moduleColors <- merge$colors
MEs <- merge$newMEs

# 准备性状数据（数值化）
# 分类变量需要编码为0/1
nSamples <- nrow(datExpr)

# 计算模块特征基因
MEs0 <- moduleEigengenes(datExpr, moduleColors)$eigengenes
MEs <- orderMEs(MEs0)

# 计算模块-性状相关性
moduleTraitCor <- cor(MEs, traitData, use = "p")
moduleTraitPvalue <- corPvalueStudent(moduleTraitCor, nSamples)

# 热图可视化
textMatrix <- paste(signif(moduleTraitCor, 2), "\n(",
                    signif(moduleTraitPvalue, 1), ")", sep = "")
dim(textMatrix) <- dim(moduleTraitCor)

par(mar = c(6, 8.5, 3, 3))
labeledHeatmap(Matrix = moduleTraitCor,
               xLabels = names(traitData),
               yLabels = names(MEs),
               ySymbols = names(MEs),
               colorLabels = FALSE,
               colors = blueWhiteRed(50),
               textMatrix = textMatrix,
               setStdMargins = FALSE,
               cex.text = 0.5,
               zlim = c(-1, 1),
               main = "Module-Trait Relationships")

# 找显著关联的模块
sig_modules <- which(moduleTraitPvalue[, "Disease"] < 0.05)
cat("Significant modules for Disease:\n")
print(names(MEs)[sig_modules])

# 计算Gene Significance和Module Membership
# 以"Disease"性状为例
trait_of_interest <- as.data.frame(traitData$Disease)
names(trait_of_interest) <- "Disease"

# GS: 基因与性状的相关性
GS <- as.data.frame(cor(datExpr, trait_of_interest, use = "p"))
names(GS) <- paste("GS.", names(trait_of_interest), sep = "")
GSPvalue <- as.data.frame(corPvalueStudent(as.matrix(GS), nSamples))

# MM: 基因与各模块ME的相关性
MM <- as.data.frame(cor(datExpr, MEs, use = "p"))
names(MM) <- paste("MM", substring(names(MEs), 3), sep = "")
MMPvalue <- as.data.frame(corPvalueStudent(as.matrix(MM), nSamples))

# 选择感兴趣的模块（如与Disease最相关的blue模块）
target_module <- "blue"

# 模块内基因
module_genes <- colnames(datExpr)[moduleColors == target_module]
cat("Genes in", target_module, "module:", length(module_genes), "\n")

# 筛选Hub基因
mm_col <- paste("MM", target_module, sep = "")
gs_col <- paste("GS.", "Disease", sep = "")

hub_criteria <- abs(MM[module_genes, mm_col]) > 0.8 &
                abs(GS[module_genes, gs_col]) > 0.2

hub_genes <- module_genes[hub_criteria]
cat("Hub genes:", length(hub_genes), "\n")
print(hub_genes)

# GS vs MM散点图
verboseScatterplot(abs(MM[module_genes, mm_col]),
                   abs(GS[module_genes, gs_col]),
                   xlab = paste("Module Membership in", target_module, "module"),
                   ylab = paste("Gene significance for Disease"),
                   main = paste("MM vs. GS\n"),
                   cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2,
                   col = target_module)

# 导出网络到Cytoscape
# 提取模块内top连接
TOM_module <- TOM[moduleColors == target_module, moduleColors == target_module]
dimnames(TOM_module) <- list(module_genes, module_genes)

# 导出top connections
cyt <- exportNetworkToCytoscape(TOM_module,
  edgeFile = paste("CytoscapeInput-edges-", target_module, ".txt", sep = ""),
  nodeFile = paste("CytoscapeInput-nodes-", target_module, ".txt", sep = ""),
  weighted = TRUE,
  threshold = 0.02,  # 连接权重阈值
  nodeNames = module_genes,
  nodeAttr = moduleColors[moduleColors == target_module])

library(clusterProfiler)
library(org.Hs.eg.db)

# 对目标模块做GO富集
module_genes_entrez <- bitr(hub_genes,
                             fromType = "SYMBOL",
                             toType = "ENTREZID",
                             OrgDb = org.Hs.eg.db)

# GO富集
ego <- enrichGO(gene = module_genes_entrez$ENTREZID,
                OrgDb = org.Hs.eg.db,
                ont = "BP",
                pAdjustMethod = "BH",
                pvalueCutoff = 0.05,
                readable = TRUE)

dotplot(ego, showCategory = 20, title = paste(target_module, "module GO enrichment"))

# KEGG富集
ekegg <- enrichKEGG(gene = module_genes_entrez$ENTREZID,
                     organism = "hsa",
                     pvalueCutoff = 0.05)

dotplot(ekegg, showCategory = 15, title = paste(target_module, "module KEGG enrichment"))

# 对所有模块批量做富集
all_modules <- unique(moduleColors)
module_enrichment <- list()

for (mod in all_modules) {
  mod_genes <- colnames(datExpr)[moduleColors == mod]
  mod_entrez <- bitr(mod_genes, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)
  if (nrow(mod_entrez) > 10) {
    module_enrichment[[mod]] <- enrichGO(gene = mod_entrez$ENTREZID,
                                          OrgDb = org.Hs.eg.db, ont = "BP",
                                          pvalueCutoff = 0.05, readable = TRUE)
  }
}

# ===== WGCNA完整流程（可直接运行）=====
library(WGCNA)
options(stringsAsFactors = FALSE)
enableWGCNAThreads(nThreads = 8)

# 1. 读取数据（行=样本, 列=基因）
datExpr <- read.csv("expr_matrix_samples_as_rows.csv", row.names = 1)
traitData <- read.csv("clinical_data.csv", row.names = 1)

# 2. 质控
gsg <- goodSamplesGenes(datExpr, verbose = 3)
if (!gsg$allOK) datExpr <- datExpr[gsg$goodSamples, gsg$goodGenes]

# 3. 软阈值
powers <- c(1:10, seq(12, 20, 2))
sft <- pickSoftThreshold(datExpr, powerVector = powers, networkType = "signed", verbose = 5)
softPower <- sft$powerEstimate
if (is.na(softPower)) softPower <- 8

# 4. 网络构建
net <- blockwiseModules(datExpr, power = softPower, networkType = "signed",
  TOMType = "signed", minModuleSize = 30, mergeCutHeight = 0.25,
  numericLabels = TRUE, saveTOMs = TRUE, verbose = 3, maxBlockSize = 20000)
moduleColors <- labels2colors(net$colors)

# 5. 模块-性状关联
MEs <- moduleEigengenes(datExpr, moduleColors)$eigengenes
MEs <- orderMEs(MEs)
moduleTraitCor <- cor(MEs, traitData, use = "p")
moduleTraitPvalue <- corPvalueStudent(moduleTraitCor, nrow(datExpr))

# 6. Hub基因
# (选择与目标性状最相关的模块)
target_module <- "blue"  # 根据热图结果选择
module_genes <- colnames(datExpr)[moduleColors == target_module]
MM <- cor(datExpr[, module_genes], MEs[, paste0("ME", target_module)], use = "p")
GS <- cor(datExpr[, module_genes], traitData$Disease, use = "p")
hub_genes <- module_genes[abs(MM) > 0.8 & abs(GS) > 0.2]