一、加载包以及文件预处理

# ====================== 1. 加载所需包 ======================
library(WGCNA)
library(tidyverse)
library(cowplot)
library(FactoMineR)
library(factoextra)
library(gplots)
library(ggpubr)
library(grid)
library(gridExtra)
library(dplyr)
library(pheatmap)
library(ggplot2)
# 全局设置
options(stringsAsFactors = FALSE)
allowWGCNAThreads() # 启用多线程加速
cat("====================== 包加载与全局设置完成 ======================\n")# ====================== 2. 设置工作目录 & 读取数据 ======================
setwd("C:/Users/27928/Desktop/长寿机制论文/转录组数据分析/分析-WGCNA") # 统一工作目录
# 读取FPKM数据（3个样本一年生叶片）
fpkm00 <- read.table("2个样本二年生枝条/2samples_2_S.txt", header = T, sep = "\t")

# 统计重复基因并提示
dup_gene <- table(duplicated(fpkm00$gene_id))
cat("重复基因统计结果：", paste(names(dup_gene), "=", dup_gene, collapse = " | "), "\n")
cat("====================== 工作目录设置与数据读取完成 ======================\n")
# ====================== 3. 数据预处理 ======================
# 将gene_id设为行名，去重（保留首次出现的基因）
fpkm <- fpkm00[!duplicated(fpkm00$gene_id), ] %>% column_to_rownames("gene_id")
# 对数转换（FPKM+1避免log(0)，防止数值错误）
data <- log2(fpkm + 1)

# 筛选基因：选择MAD最高的前5000个基因（可根据需求修改为10000/15000）
keep_num <- 5000 
keep_data <- data[order(apply(data, 1, mad), decreasing = T)[1:keep_num], ]

# 转置数据（WGCNA强制要求：行=样本，列=基因）
datExpr0 <- as.data.frame(t(keep_data))
nGenes <- ncol(datExpr0)
nSamples <- nrow(datExpr0)
cat("预处理后：样本数 =", nSamples, "，筛选后基因数 =", nGenes, "\n")

# 保存预处理后的表达矩阵（用于后续复用）
save(datExpr0, file = "fpkm_trait-dataInput.RData") # 与模块可视化代码的加载文件对应
cat("====================== 数据预处理完成，已保存预处理表达矩阵 ======================\n")

二、进行数据分析

# 核心：匹配样本名，设置Y_L(幼龄)、A_1_L(成年)、O_1_L(老龄)分组
sample_names <- colnames(data)
datTraits <- data.frame(
  row.names = sample_names,
  group = c("A_2_S_1","A_2_S_2","A_2_S_3","O_2_S_1","O_2_S_2","O_2_S_3")
)
# 给分组添加数字编号（用于可视化上色，避免字符型无法绘图）
grouptype <- data.frame(group = sort(unique(datTraits$group)), groupNo = 1:length(unique(datTraits$group)))
datTraits$groupNo <- "NA"
for(i in 1:nrow(grouptype)){
  datTraits[which(datTraits$group == grouptype$group[i]), 'groupNo'] <- grouptype$groupNo[i]
}
datTraits$groupNo <- as.numeric(datTraits$groupNo) # 转为数值型
# 保存性状数据
write.table(datExpr0, file = "datExpr0.txt", sep = "\t", quote = F)
write.csv(datTraits, file = "datTraits.txt", row.names = T, quote = F)
cat("====================== 性状数据构建完成，已保存样本分组信息 ======================\n")
# ====================== 5. 数据质量控制 ======================
# 5.1 缺失值 & 低质量基因/样本过滤（WGCNA核心质控步骤）
gsg <- goodSamplesGenes(datExpr0, verbose = 3)
if (!gsg$allOK) {
  # 过滤低质量基因
  if (sum(!gsg$goodGenes) > 0) {
    printFlush(paste("移除低质量基因数量：", sum(!gsg$goodGenes)))
  }
  # 过滤离群样本并同步更新性状数据
  if (sum(!gsg$goodSamples) > 0) {
    printFlush(paste("移除离群样本数量：", sum(!gsg$goodSamples)))
    datTraits <- datTraits[gsg$goodSamples, ] 
  }
  datExpr0 <- datExpr0[gsg$goodSamples, gsg$goodGenes]
}
cat("质控后：样本数 =", nrow(datExpr0), "，基因数 =", ncol(datExpr0), "\n")

# 5.2 样本聚类树（核心修正：按Y_L/A_1_L/O_1_L三组上色，同组同色）- 可视化样本分群情况
sampleTree <- hclust(dist(datExpr0), method = "average") # 保留原聚类计算，无需修改

# 核心步骤1：重构分组逻辑，提取每个样本的「组前缀」（Y_L/A_1_L/O_1_L）
# 从datTraits$group中提取前3/5个字符，得到三组标签（适配Y_L/A_1_L/O_1_L命名）
datTraits$group_3class <- substr(datTraits$group, 1, ifelse(substr(datTraits$group, 1, 2) == "A_", 3, 5))
# 手动指定三组顺序（幼龄→成年→老龄），避免字母排序混乱
datTraits$group_3class <- factor(datTraits$group_3class, levels = c("A_2_S", "O_2_S"))
# 给三组分配数字编号（用于颜色映射，1=Y_L,2=A_1_L,3=O_1_L）
datTraits$group_3class_no <- as.numeric(datTraits$group_3class)

# 核心步骤2：定义三组专属配色（论文级美观配色，同组同色，对比度高）
# 配色方案：Y_L(幼龄)=蓝色，A_1_L(成年)=橙色，O_1_L(老龄)=红色（可根据需求修改）
three_class_colors <- c("#4E79A7", "#E15759") 
# 按三组编号映射颜色，实现9个样本→3种颜色（同组同色）
sample_colors <- three_class_colors[datTraits$group_3class_no]

# 绘图并保存（优化尺寸+边距，适配三组标注，防止标签重叠）
graphics.off() # 关闭多余绘图设备，防止报错
pdf("step1_Sample_dendrogram_trait_3class.pdf", width = 14, height = 10) # 加宽宽度，适配样本名显示
par(mar = c(2, 4, 3, 1), cex = 0.9) # 调整下边缘距，避免样本名被遮挡，放大字体
plotDendroAndColors(
  sampleTree, 
  sample_colors, 
  #groupLabels = "Y_L/A_1_L/O_1_L Group", # 图例标题明确标注“三组”
  cex.dendroLabels = 0.9, # 放大样本名标签，提升可读性
  marAll = c(2, 4, 3, 1), # 统一调整所有子图边距，匹配par设置
  main = "Sample Dendrogram (A_1_L:Adult, O_1_L:Old)", # 标题标注分组含义
  cex.main = 1.2 # 放大主标题
)
dev.off()

# 打印分组提示，验证分组结果（运行时可查看，确保无错误）
cat("样本三组分组完成：\n")
print(table(datTraits$group, datTraits$group_3class)) # 显示每个样本对应的三组标签
cat("5.2 样本聚类树绘制完成：按Y_L/A_1_L/O_1_L三组上色，结果保存为step1_Sample_dendrogram_trait_3class.pdf\n")
# 5.3 PCA分析验证样本分组（辅助质控，验证分组合理性）
group_list <- datTraits$groupNo
dat.pca <- PCA(datExpr0, graph = F)
pca <- fviz_pca_ind(
  dat.pca,
  title = "Principal Component Analysis (PCA)",
  legend.title = "Groups",
  geom.ind = c("point","text"),
  pointsize = 2,
  labelsize = 4,
  repel = TRUE, # 避免标签重叠
  col.ind = group_list,
  axes.linetype = NA,
  mean.point = F
) + theme(legend.position = "none") + coord_fixed(ratio = 1)
ggsave("step1_PCA.pdf", pca, width = 8, height = 6)
print(pca)
cat("====================== 数据质量控制完成，已生成样本聚类树和PCA图 ======================\n")

# ====================== 6. 选择最佳软阈值（构建无标度网络） ======================
datExpr <- datExpr0 # 统一表达矩阵变量名
R.sq_cutoff <- 0.85 # 无标度网络拟合度阈值（常规阈值0.85）
powers <- c(seq(1,20,by = 1), seq(22,30,by = 2)) # 软阈值候选范围

# 计算软阈值（耗时步骤，多线程已加速）
sft <- pickSoftThreshold(
  datExpr,
  networkType = "unsigned", # 无符号网络（仅关注相关性大小，不关注正负）
  powerVector = powers,
  RsquaredCut = R.sq_cutoff,
  verbose = 5
)

# 绘制软阈值选择图（保存为PDF，便于查看）
graphics.off()
pdf("step2_power-value.pdf", width = 16, height = 12)
par(mfrow = c(1,2)); cex1 = 0.9
# 左图：软阈值与无标度拟合度的关系
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")
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     labels = powers, cex = cex1, col = "red")
abline(h = R.sq_cutoff, col = "red", lwd = 2, lty = 2) # 拟合度阈值线
# 右图：软阈值与平均连接度的关系
plot(sft$fitIndices[,1], sft$fitIndices[,5],
     xlab = "Soft Threshold (power)", ylab = "Mean Connectivity", type = "n")
text(sft$fitIndices[,1], sft$fitIndices[,5], labels = powers, cex = cex1, col = "red")
abline(h = 100, col = "red", lwd = 2, lty = 2) # 平均连接度参考线
dev.off()

# 确定最终power值（针对小样本优化，无自动估计时设默认值12）
power <- sft$powerEstimate
if(is.na(power)){
  power <- 12  # 9个样本推荐默认值12，适配小样本网络构建
}
cat("最终选择的软阈值power =", power, "\n")
cat("====================== 软阈值选择完成，已生成软阈值分析图 ======================\n")

# ====================== 7. 构建WGCNA共表达网络 & 识别模块 ======================
net <- blockwiseModules(
  datExpr,
  power = power, # 上述选择的软阈值
  maxBlockSize = 5000, # 单块最大基因数，匹配筛选的5000基因
  corType = "pearson", # 皮尔森相关系数，常规选择
  networkType = "unsigned",
  TOMType = "unsigned", # TOM矩阵类型与网络类型一致
  minModuleSize = 20,  # 小样本降低最小模块大小，保证模块数量
  mergeCutHeight = 0.25, # 模块合并阈值（高度<0.25则合并，值越小模块越多）
  numericLabels = TRUE, # 先以数字标记模块，后续转颜色
  saveTOMs = TRUE,  # 开启TOM矩阵保存（用于后续网络可视化导出）
  saveTOMFileBase = "WGCNA_TOM",  # TOM矩阵保存前缀
  verbose = 3
)

# 保存网络构建结果（与模块可视化代码的加载文件对应）
save(net, datExpr, datTraits, power, file = "networkConstruction-stepByStep.RData")

# 模块数字编号转颜色标签（WGCNA标准颜色命名，便于可视化）
moduleColors <- labels2colors(net$colors)
module_count <- table(moduleColors)
cat("模块数量及各模块基因数：\n")
print(module_count)

# 绘制基因聚类树+模块颜色（核心可视化，展示基因分群结果）
graphics.off()
pdf("step3_gene_dendrogram_module_colors.pdf", width = 16, height = 12)
plotDendroAndColors(
  net$dendrograms[[1]], 
  moduleColors[net$blockGenes[[1]]],
  "Module colors",
  dendroLabels = FALSE, # 关闭基因名标签，避免图面杂乱
  hang = 0.03,
  addGuide = TRUE, 
  guideHang = 0.05,
  main = "Gene dendrogram and module colors"
)
dev.off()
cat("====================== WGCNA网络构建完成，已生成基因聚类树+模块颜色图 ======================\n")

# ====================== 8. 模块与样本分组性状的关联分析 ======================
# 8.1 计算模块特征基因（MEs）：每个模块的代表性基因，反映模块整体表达趋势
MEs <- moduleEigengenes(datExpr, moduleColors)$eigengenes
# 移除grey模块（灰色模块为未分群基因，无生物学意义）
if ("MEgrey" %in% colnames(MEs)) {
  MEs <- MEs[, -which(colnames(MEs) == "MEgrey")]
}
# 异常检查：确保有有效模块生成
if(ncol(MEs) == 0){
  stop("错误：无有效模块生成！请降低minModuleSize或调整软阈值power！")
}

# 8.2 分组性状转数值矩阵（哑变量，用于相关性计算）
trait_matrix <- model.matrix(~ 0 + datTraits$group) %>% 
  `colnames<-`(levels(factor(datTraits$group))) %>% 
  as.data.frame() %>% 
  `rownames<-`(rownames(datTraits))

# 8.3 计算模块-性状相关性 & 显著性P值
moduleTraitCor <- cor(MEs, trait_matrix, use = "p") # 皮尔森相关
moduleTraitPvalue <- corPvalueStudent(moduleTraitCor, nrow(datExpr)) # 学生t检验计算P值

# 8.4 可视化：模块-性状相关性热图（指定性状顺序，关闭自动聚类）
# 构建带P值的文本矩阵（用于热图标注）
text_matrix <- paste(
  signif(moduleTraitCor, 2), 
  "\n(", 
  signif(moduleTraitPvalue, 1), 
  ")", 
  sep = ""
)
dim(text_matrix) <- dim(moduleTraitCor)
dimnames(text_matrix) <- dimnames(moduleTraitCor) # 关键：添加维度名，避免绘图报错

# 构建行注释（模块颜色，便于识别）
row_anno <- data.frame(
  ModuleColor = gsub("ME", "", rownames(moduleTraitCor)),  # 提取模块颜色名
  row.names = rownames(moduleTraitCor)
)
anno_colors <- list(
  ModuleColor = setNames(unique(gsub("ME", "", rownames(moduleTraitCor))), 
                         unique(gsub("ME", "", rownames(moduleTraitCor))))
)

# 按指定顺序排序（模块按系统发育树，性状按Y/A/O分组）
# 1. 模块排序：基于MEs聚类的系统发育树顺序
module_dist <- dist(t(MEs), method = "euclidean")
module_tree <- hclust(module_dist, method = "average")
module_order <- module_tree$labels[module_tree$order]
# 2. 性状排序：指定幼龄→成年→老龄的固定顺序
trait_order <- c("A_2_S_1","A_2_S_2","A_2_S_3","O_2_S_1","O_2_S_2","O_2_S_3")
trait_order <- trait_order[trait_order %in% colnames(moduleTraitCor)] # 安全过滤，匹配实际列名

# 调整矩阵顺序（匹配指定排序）
moduleTraitCor_ordered <- moduleTraitCor[module_order, trait_order, drop = FALSE]
text_matrix_ordered <- text_matrix[module_order, trait_order, drop = FALSE]
row_anno_ordered <- row_anno[module_order, , drop = FALSE]

# 绘制并保存排序后的相关性热图
pheatmap(
  moduleTraitCor_ordered,
  color = colorRampPalette(c("blue", "white", "red"))(100), # 蓝-白-红配色，负相关→正相关
  display_numbers = text_matrix_ordered,
  number_color = "black",
  fontsize_number = 5,
  fontsize_row = 8,
  fontsize_col = 8,
  cluster_rows = FALSE,  # 关闭行聚类，使用指定的系统发育树顺序
  cluster_cols = FALSE,  # 关闭列聚类，使用指定的样本分组顺序
  main = "Module-Trait Relationships (Y→A→O Order)",
  annotation_row = row_anno_ordered,
  annotation_colors = anno_colors,
  show_rownames = TRUE,
  show_colnames = TRUE,
  filename = "step4_Module-Trait_Correlation_Specified_Order.pdf",
  width = 12,
  height = 8
)

# 8.5 筛选关键模块（按各性状的相关性绝对值最大筛选）
key_modules <- apply(moduleTraitCor_ordered, 2, function(x) {
  module_idx <- which.max(abs(x))
  module_name <- rownames(moduleTraitCor_ordered)[module_idx]
  gsub("ME", "", module_name)
})
cat("关键模块与样本分组的对应关系（指定性状顺序）：\n")
print(key_modules)
cat("====================== 模块-性状关联分析完成，已生成相关性热图 ======================\n")

# ====================== 9. 结果保存（核心结果文件，便于后续分析） ======================
# 模块-基因映射表（所有基因的模块注释）
module_info <- data.frame(
  Gene = colnames(datExpr),
  Module = moduleColors,
  row.names = colnames(datExpr),
  check.names = FALSE
)
# 保存各类结果文件
write.csv(module_info, "module_gene_mapping.csv", row.names = TRUE, quote = FALSE)
write.csv(as.data.frame(moduleTraitCor_ordered), "module_trait_correlation_specified_order.csv", row.names = TRUE, quote = FALSE)
write.csv(as.data.frame(moduleTraitPvalue), "module_trait_pvalue.csv", row.names = TRUE, quote = FALSE)
cat("已保存核心结果文件：\n")
cat("  - 模块-基因映射表：module_gene_mapping.csv\n")
cat("  - 模块-性状相关性矩阵：module_trait_correlation_specified_order.csv\n")
cat("  - 模块-性状相关性P值：module_trait_pvalue.csv\n")
cat("====================== 核心结果保存完成 ======================\n")

# ====================== 10. 模块基因表达热图 + 特征基因（ME）表达量可视化【固定顺序+热图无样本名+柱状图样本名防重叠】 ======================
####### 不同模块基因表达热图及特征基因（ME）表达量分析（固定样本顺序+热图无样本名+柱状图显样本名无重叠）
cat("\n===== 开始绘制各模块基因表达热图 + 特征基因表达量（固定顺序+热图无样本名） =====\n")
# 加载网络数据（与原代码加载逻辑一致，确保变量匹配）
lnames = load(file = "networkConstruction-stepByStep.RData");
# 计算基因表达相关性矩阵（皮尔森）和TOM矩阵（用于验证，可选）
person_cor <- cor(datExpr0, use = 'p') # 基因-基因表达相关性
TOM_corr <- TOMsimilarityFromExpr(datExpr, power = power, TOMType = "unsigned") # 重构TOM矩阵
Colors <- moduleColors # 模块颜色标签
# 统一矩阵行/列名和颜色标签命名（避免匹配错误）
colnames(TOM_corr) <- rownames(TOM_corr) <- colnames(datExpr0)
colnames(person_cor) <- rownames(person_cor) <- colnames(datExpr0)
names(Colors) <- colnames(datExpr0)

# 创建结果文件夹（存放模块可视化图，避免文件杂乱）
out_dir <- "15_module_displaying module heatmap and the eigengene expression"
dir.create(out_dir, showWarnings = FALSE, recursive = TRUE)

# 提取有效模块（排除grey模块，无生物学意义）
umc <- unique(Colors)[unique(Colors) != "grey"]
lumc <- length(umc)
cat("有效模块数量（排除grey）：", lumc, "，模块列表：", paste(umc, collapse = ", "), "\n")

# 核心步骤1：手动定义固定样本顺序（严格按需求）
sample_fixed_order <- c("A_2_S_1","A_2_S_2","A_2_S_3","O_2_S_1","O_2_S_2","O_2_S_3")
# 安全过滤：仅保留数据中实际存在的样本名（避免报错）
sample_fixed_order <- sample_fixed_order[sample_fixed_order %in% rownames(datExpr0)]
cat("最终绘图样本固定顺序：", paste(sample_fixed_order, collapse = " → "), "\n")

# 核心步骤2：按固定顺序重新排列表达矩阵和ME矩阵
datExpr0_ordered <- datExpr0[sample_fixed_order, ]
MEs_ordered <- MEs[sample_fixed_order, ]

# 循环绘制每个模块的「基因表达热图（无样本名） + ME表达量柱状图（显样本名无重叠）」
for (i in 1:lumc) {
  module <- umc[i]
  ME <- MEs_ordered[, paste("ME", module, sep = "")] # 提取排序后的模块ME值
  # 定义输出PDF文件名
  pdf_file <- paste0(out_dir, "/15_", module, "_module_heatmap_eigengene_expression.pdf")
  # 绘制双图：加宽图幅适配柱状图样本名，高度微调增加边距空间
  pdf(file = pdf_file, width = 10, height = 8) 
  # 上图：模块基因表达热图（固定顺序+无样本名，严格按需求）
  par(mfrow = c(2, 1), mar = c(0.3, 5.5, 3, 2)) # 保留原边距，无样本名无需调整
  plotMat(
    t(scale(datExpr0_ordered[, Colors == module])), # 行标准化，便于跨基因比较
    nrgcols = 30, # 颜色梯度数量，渐变平滑
    rlabels = FALSE, # 关闭行标签（基因名）
    clabels = FALSE, # 核心：关闭列标签（样本名），严格遵循需求
    rcols = module, # 热图颜色匹配模块颜色
    main = paste(module, "Module (", sum(Colors == module), "Genes)", sep = " "), # 标题含模块基因数
    cex.main = 2
  )
  # 下图：特征基因（ME）表达量柱状图（显样本名+解决与xlab重叠问题）
  par(mar = c(8, 4.2, 0, 0.7)) # 核心修正：将下边缘距从5增至8，为样本名预留足够空间
  barplot(
    ME, 
    col = module, # 柱状图颜色匹配模块颜色
    names.arg = sample_fixed_order, # 横轴添加固定顺序的样本名标签
    las = 2, # 样本名垂直倾斜90°，防止横向重叠
    cex.names = 0.8, # 样本名字体大小
    ylab = "Eigengene Expression (ME)", # 纵轴标签：特征基因表达量
    #xlab = "Sample", # 横轴标签
    cex.axis = 0.9, # 坐标轴刻度字体大小
    mgp = c(3, 1, 0) # 优化坐标轴标签位置，进一步避免重叠
  )
  dev.off()
  cat("  已生成", module, "模块可视化图（热图无样本名+柱状图样本名无重叠）：", pdf_file, "\n")
}

cat("===== 所有模块热图+ME表达量可视化完成，结果存放于：", out_dir, "=====\n")
cat("====================== 模块可视化分析完成 ======================\n")

# ====================== 11. 导出模块基因网络文件（VisANT/Cytoscape） ======================
####### gene export for VisANT or cytoscape（网络可视化文件导出）
if(T){  
  # 可自定义参数（根据关键模块调整）
  target_module = "brown"  # 要导出的模块（可改为key_modules中的关键模块，如"blue"/"turquoise"）
  top_gene_num = 100       # 筛选Top N核心基因（按软连接度排序）
  cytoscape_threshold = 0.02  # Cytoscape边过滤阈值（>0.02保留，过滤低相关性边）
  
  # 加载TOM矩阵并命名
  TOM <- TOMsimilarityFromExpr(datExpr, power = power, TOMType = "unsigned")
  dimnames(TOM) <- list(colnames(datExpr), colnames(datExpr))
  
  # 异常检查：确保关键变量存在
  if(!exists("datExpr")) stop("错误：未找到表达矩阵datExpr！")
  if(!exists("moduleColors")) stop("错误：未找到模块注释moduleColors！")
  if(!exists("TOM")) stop("错误：未找到TOM矩阵，请确认saveTOMs=TRUE！")
  
  # 提取目标模块基因
  probes = colnames(datExpr)
  inModule = (moduleColors == target_module)
  modProbes = probes[inModule]
  if(length(modProbes) == 0){stop(paste("错误：模块", target_module, "无基因！"))}
  
  # 筛选Top N核心基因（按软连接度）
  modTOM = TOM[inModule, inModule]
  dimnames(modTOM) = list(modProbes, modProbes)
  IMConn = softConnectivity(datExpr[, modProbes])
  top = (rank(-IMConn) <= top_gene_num)
  filterTOM = modTOM[top, top]
  topProbes = modProbes[top]
  
  # 导出VisANT格式文件
  vis_file = paste("visANTinput-", target_module, ".txt", sep = "")
  exportNetworkToVisANT(filterTOM, file = vis_file, weighted = T, threshold = 0)
  # 导出Cytoscape格式文件（边+节点）
  edge_file = paste("CytoscapeInput-edges-", target_module, ".txt", sep = "")
  node_file = paste("CytoscapeInput-nodes-", target_module, ".txt", sep = "")
  exportNetworkToCytoscape(filterTOM, edgeFile = edge_file, nodeFile = node_file,
                           weighted = TRUE, threshold = cytoscape_threshold,
                           nodeNames = topProbes, nodeAttr = moduleColors[inModule][top])
  
  cat("\n===== 网络可视化文件导出完成 =====\n")
  cat("  - VisANT文件：", vis_file, "\n")
  cat("  - Cytoscape边文件：", edge_file, "\n")
  cat("  - Cytoscape节点文件：", node_file, "\n")
}

# ====================== 所有分析流程结束 ======================
cat("\n====================== 全部WGCNA分析流程执行完成！ ======================\n")
cat("生成的结果文件分类：\n")
cat("1. 质控与软阈值分析：step1_*.pdf、step2_power-value.pdf\n")
cat("2. 网络构建与模块分群：step3_gene_dendrogram_module_colors.pdf\n")
cat("3. 模块-性状关联：step4_Module-Trait_Correlation_Specified_Order.pdf、*.csv\n")
cat("4. 模块可视化：", out_dir, "/（各模块热图+ME表达量PDF）\n")
cat("5. 网络可视化文件：visANTinput-*.txt、CytoscapeInput-*.txt\n")

模块与表达量的结果文件