Bulk RNA-seq，cpm表达矩阵和分组信息被放在2个不同的txt中，使用WGCNA包分析基因表达和性状的关联性，并可视化结果，R语言分析代码

好的，以下是基于WGCNA包的Bulk RNA-seq数据分析和可视化的R语言代码：

# 导入所需的包
library(WGCNA)

# 读取cpm表达矩阵和分组信息
expr_matrix <- as.data.frame(read.table("cpm_matrix.txt", header = TRUE, row.names = 1))
group_info <- read.table("group_info.txt", header = TRUE)

# 将分组信息转换为因子型
group_info$group <- as.factor(group_info$group)

# 构建WGCNA网络
powers <- c(1:10)
sft <- pickSoftThreshold(expr_matrix, powerVector = powers, verbose = 5)
power <- sft$powerEstimate
net <- blockwiseModules(expr_matrix, power = power, TOMType = "unsigned",
                        minModuleSize = 30, mergeCutHeight = 0.25,
                        numericLabels = TRUE, pamRespectsDendro = FALSE,
                        saveTOMs = TRUE, saveTOMFileBase = "TOM", verbose = 0)

# 可视化模块-性状关联
module_trait <- moduleTraitCor(net, traitData = group_info, usePearson = TRUE, verbose = 0)
module_colors <- labels2colors(net$colors)
plotModuleTrait(net, moduleTrait, labelModule = module_colors, xaxisLabel = "Module Membership", yaxisLabel = "Trait Correlation")

# 可视化模块-基因关联
module_genes <- moduleEigengenes(expr_matrix, colors = net$colors, verbose = 0)
ME_list <- orderMEs(module_genes)
for (k in 1:length(ME_list)) {
  ME <- ME_list[k]
  module_color <- module_colors[ME]
  module_genes <- module_genes[ME, ]
  module_genes <- module_genes[order(module_genes, decreasing = TRUE), , drop = FALSE]
  top_genes <- names(module_genes)[1:10]
  pdf(paste("ME", ME, "-", module_color, ".pdf", sep = ""))
  barplot(module_genes, main = paste("Module Eigengene Expression Profile in", module_color, "Module", sep = "\n"))
  dev.off()
}

这里，我们读取了cpm表达矩阵和分组信息，并将分组信息转换为因子型。然后，我们使用`pickSoftThreshold()`函数选择最佳的WGCNA网络构建参数，使用`blockwiseModules()`函数构建WGCNA网络并保存TOM矩阵。接着，我们使用`moduleTraitCor()`函数可视化模块-性状关联，并使用`plotModuleTrait()`函数绘制模块-性状关联图。最后，我们使用`moduleEigengenes()`函数计算模块的基因表达量，并使用`barplot()`函数绘制模块-基因关联图。

请根据实际情况修改代码中的文件路径和参数，并根据需要调整可视化图表的标题和标签。