Trait: drought severity
        |
        | corr = 0.82, p = 1e-05
        v
turquoise module eigengene
        |
        | high module membership
        v
Hub genes: NAC1, DREB2A, HSP70, ...

flowchart TD
    A[RNA-seq count / TPM matrix] --> B[VST / log2(TPM+1) 转换]
    B --> C[低表达和低变异基因过滤]
    C --> D[样本聚类检测离群样本]
    D --> E[选择 soft-threshold power]
    E --> F[构建 adjacency matrix]
    F --> G[计算 TOM similarity]
    G --> H[层次聚类识别 modules]
    H --> I[合并相似 modules]
    I --> J[计算 module eigengenes]
    J --> K[模块-性状相关热图]
    K --> L[筛选 hub genes]
    L --> M[GO/KEGG 富集和 Cytoscape 网络]

wgcna_project/
├── 00_metadata/
│   └── traits.csv
├── 01_expression/
│   ├── raw_counts.csv
│   ├── vst_matrix.csv
│   └── filtered_expr.csv
├── 02_qc/
│   ├── sample_clustering.pdf
│   └── soft_power.pdf
├── 03_wgcna/
│   ├── network.RData
│   ├── module_genes.csv
│   ├── module_trait_correlation.csv
│   └── hub_genes.csv
├── 04_enrichment/
├── 05_cytoscape/
└── report/

行 = 样本
列 = 基因

library(DESeq2)
library(tidyverse)

count_data <- read.csv("01_expression/raw_counts.csv", row.names = 1, check.names = FALSE)
trait_data <- read.csv("00_metadata/traits.csv", row.names = 1)

count_data <- count_data[, rownames(trait_data)]

dds <- DESeqDataSetFromMatrix(
  countData = round(as.matrix(count_data)),
  colData = trait_data,
  design = ~ 1
)

keep <- rowSums(counts(dds) >= 10) >= 3
dds <- dds[keep, ]

vsd <- vst(dds, blind = TRUE)
vst_matrix <- assay(vsd)

write.csv(vst_matrix, "01_expression/vst_matrix.csv")

tpm <- read.csv("01_expression/tpm_matrix.csv", row.names = 1, check.names = FALSE)
log_tpm <- log2(tpm + 1)
write.csv(log_tpm, "01_expression/log2_tpm_matrix.csv")

expr <- read.csv("01_expression/vst_matrix.csv", row.names = 1, check.names = FALSE)

gene_mad <- apply(expr, 1, mad)
top_genes <- names(sort(gene_mad, decreasing = TRUE))[1:min(8000, length(gene_mad))]
expr_filtered <- expr[top_genes, ]

write.csv(expr_filtered, "01_expression/filtered_expr.csv")

library(WGCNA)

options(stringsAsFactors = FALSE)
allowWGCNAThreads()

expr <- read.csv("01_expression/filtered_expr.csv", row.names = 1, check.names = FALSE)
datExpr <- as.data.frame(t(expr))

gsg <- goodSamplesGenes(datExpr, verbose = 3)
if (!gsg$allOK) {
  datExpr <- datExpr[gsg$goodSamples, gsg$goodGenes]
}

sample_tree <- hclust(dist(datExpr), method = "average")

pdf("02_qc/sample_clustering.pdf", width = 10, height = 6)
plot(sample_tree, main = "Sample clustering to detect outliers", sub = "", xlab = "")
dev.off()

powers <- c(1:10, seq(12, 30, by = 2))

sft <- pickSoftThreshold(
  datExpr,
  powerVector = powers,
  networkType = "signed",
  verbose = 5
)

pdf("02_qc/soft_power.pdf", width = 10, height = 5)
par(mfrow = c(1, 2))

plot(
  sft$fitIndices[, 1],
  -sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2],
  xlab = "Soft Threshold power",
  ylab = "Scale Free Topology Model Fit, signed R^2",
  type = "n",
  main = "Scale independence"
)
text(
  sft$fitIndices[, 1],
  -sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2],
  labels = powers,
  col = "red"
)
abline(h = 0.8, col = "red")

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, col = "red")
dev.off()

如果 power = 12 时 R² 达到 0.82，mean connectivity 仍然不低，
可以选择 softPower = 12。
如果 R² 一直达不到 0.8，选择曲线趋于平稳处，并在报告中说明。

softPower <- 12

net <- blockwiseModules(
  datExpr,
  power = softPower,
  networkType = "signed",
  TOMType = "signed",
  minModuleSize = 30,
  reassignThreshold = 0,
  mergeCutHeight = 0.25,
  numericLabels = TRUE,
  pamRespectsDendro = FALSE,
  saveTOMs = TRUE,
  saveTOMFileBase = "03_wgcna/TOM",
  verbose = 3
)

moduleColors <- labels2colors(net$colors)
MEs <- net$MEs

save(datExpr, net, moduleColors, MEs, file = "03_wgcna/network.RData")

pdf("03_wgcna/gene_dendrogram_modules.pdf", width = 12, height = 6)
plotDendroAndColors(
  net$dendrograms[[1]],
  moduleColors[net$blockGenes[[1]]],
  "Module colors",
  dendroLabels = FALSE,
  hang = 0.03,
  addGuide = TRUE,
  guideHang = 0.05
)
dev.off()

上方树状图：基因按照共表达相似性聚类。
下方颜色条：WGCNA 分配的模块颜色。
同一颜色中的基因表达模式更相似。
grey 模块通常表示未归入明确模块的基因。

trait_data <- read.csv("00_metadata/traits.csv", row.names = 1)
trait_data <- trait_data[rownames(datExpr), ]

trait_data$condition_binary <- ifelse(trait_data$condition == "Treat", 1, 0)

numeric_traits <- trait_data |>
  dplyr::select(where(is.numeric))

MEs0 <- moduleEigengenes(datExpr, moduleColors)$eigengenes
MEs <- orderMEs(MEs0)

module_trait_cor <- cor(MEs, numeric_traits, use = "p")
module_trait_pvalue <- corPvalueStudent(module_trait_cor, nSamples = nrow(datExpr))

text_matrix <- paste(
  signif(module_trait_cor, 2),
  "
(",
  signif(module_trait_pvalue, 1),
  ")",
  sep = ""
)
dim(text_matrix) <- dim(module_trait_cor)

pdf("03_wgcna/module_trait_heatmap.pdf", width = 9, height = 8)
labeledHeatmap(
  Matrix = module_trait_cor,
  xLabels = names(numeric_traits),
  yLabels = names(MEs),
  ySymbols = names(MEs),
  colorLabels = FALSE,
  colors = blueWhiteRed(50),
  textMatrix = text_matrix,
  setStdMargins = FALSE,
  cex.text = 0.8,
  zlim = c(-1, 1),
  main = "Module-trait relationships"
)
dev.off()

write.csv(module_trait_cor, "03_wgcna/module_trait_correlation.csv")
write.csv(module_trait_pvalue, "03_wgcna/module_trait_pvalue.csv")

target_module <- "turquoise"
target_trait <- "drought_score"

module_genes <- moduleColors == target_module

MM <- as.data.frame(cor(datExpr, MEs, use = "p"))
MMPvalue <- as.data.frame(corPvalueStudent(as.matrix(MM), nrow(datExpr)))

GS <- as.data.frame(cor(datExpr, numeric_traits[[target_trait]], use = "p"))
GSPvalue <- as.data.frame(corPvalueStudent(as.matrix(GS), nrow(datExpr)))

gene_info <- data.frame(
  gene_id = colnames(datExpr),
  module = moduleColors,
  module_membership = MM[[paste0("ME", target_module)]],
  module_membership_pvalue = MMPvalue[[paste0("ME", target_module)]],
  gene_significance = GS[, 1],
  gene_significance_pvalue = GSPvalue[, 1]
)

hub_genes <- gene_info |>
  dplyr::filter(
    module == target_module,
    abs(module_membership) >= 0.8,
    abs(gene_significance) >= 0.5
  ) |>
  dplyr::arrange(desc(abs(module_membership)), desc(abs(gene_significance)))

write.csv(gene_info, "03_wgcna/all_gene_module_membership.csv", row.names = FALSE)
write.csv(hub_genes, "03_wgcna/hub_genes_turquoise.csv", row.names = FALSE)

library(ggplot2)
library(ggrepel)

plot_df <- gene_info |>
  dplyr::filter(module == target_module) |>
  dplyr::mutate(label = ifelse(row_number() <= 10, gene_id, NA))

ggplot(plot_df, aes(abs(module_membership), abs(gene_significance))) +
  geom_point(alpha = 0.7, color = target_module) +
  geom_text_repel(aes(label = label), max.overlaps = 20) +
  theme_bw() +
  labs(
    x = paste0("Module Membership in ", target_module),
    y = paste0("Gene Significance for ", target_trait),
    title = "Hub gene screening"
  )

ggsave("03_wgcna/MM_vs_GS_turquoise.pdf", width = 7, height = 6)

x 轴越靠右：基因越能代表该模块。
y 轴越靠上：基因越和目标性状相关。
右上角基因：优先候选 hub genes。

load("03_wgcna/network.RData")

target_module <- "turquoise"
module_genes <- moduleColors == target_module
module_expr <- datExpr[, module_genes]

TOM <- TOMsimilarityFromExpr(
  module_expr,
  power = softPower,
  networkType = "signed",
  TOMType = "signed"
)

dimnames(TOM) <- list(colnames(module_expr), colnames(module_expr))

exportNetworkToCytoscape(
  TOM,
  edgeFile = "05_cytoscape/turquoise_edges.txt",
  nodeFile = "05_cytoscape/turquoise_nodes.txt",
  weighted = TRUE,
  threshold = 0.08,
  nodeNames = colnames(module_expr),
  nodeAttr = target_module
)

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

target_genes <- gene_info |>
  dplyr::filter(module == "turquoise") |>
  dplyr::pull(gene_id)

entrez <- mapIds(
  org.Hs.eg.db,
  keys = target_genes,
  column = "ENTREZID",
  keytype = "ENSEMBL",
  multiVals = "first"
) |>
  na.omit()

ego <- enrichGO(
  gene = entrez,
  OrgDb = org.Hs.eg.db,
  keyType = "ENTREZID",
  ont = "BP",
  pAdjustMethod = "BH",
  pvalueCutoff = 0.05
)

write.csv(as.data.frame(ego), "04_enrichment/turquoise_GO_BP.csv", row.names = FALSE)

turquoise 模块与干旱评分显著正相关，说明该模块基因整体随干旱程度升高而上调。
该模块富集于水分胁迫和 ABA 响应通路，其中 DREB2A 同时具有较高 module membership 和 gene significance，
因此可作为后续验证的候选调控基因。