372_微生物组宿主转录组联合¶

library(mixOmics)   # 多组学整合分析包
library(compositions)  # CLR转换

# 加载数据
microbiome <- read.csv("otu_relative.csv", row.names = 1)  # 样本×OTU
transcriptome <- read.csv("gene_expression.csv", row.names = 1)  # 样本×基因

# 步骤1：菌群数据 CLR 转换（消除成分数据约束）
microbiome_clr <- as.data.frame(clr(microbiome + 0.5))  # +0.5 处理零值

# 步骤2：转录组数据标准化
transcriptome_scaled <- scale(transcriptome)   # Z-score标准化

# 步骤3：sPLS 联合降维
spls_result <- spls(
  X = microbiome_clr,        # 自变量：菌群矩阵（样本×OTU）
  Y = transcriptome_scaled,  # 因变量：转录组矩阵（样本×基因）
  ncomp = 3,                 # 提取3个潜在成分
  keepX = c(10, 10, 10),     # 每个成分保留X侧10个特征（菌）
  keepY = c(20, 20, 20)      # 每个成分保留Y侧20个特征（基因）
)

# 查看关联图
plotVar(spls_result,
  comp = 1,           # 查看第1个潜在成分
  var.names = TRUE,   # 显示变量名
  cex = c(0.8, 0.6)  # 调整字体大小
)

# 提取关键特征对
selected_micro <- selectVar(spls_result, comp = 1)$X$name  # 选出的关键菌
selected_gene  <- selectVar(spls_result, comp = 1)$Y$name  # 选出的关键基因
cat("关键菌：", selected_micro, "\n")
cat("关键基因：", selected_gene, "\n")

library(WGCNA)     # 加权基因共表达网络分析
library(vegan)     # 用于 CLR 转换

# 步骤1：WGCNA 构建基因模块
enableWGCNAThreads(4)   # 开启多线程

# 自动软阈值选择
powers <- c(1:20)
sft <- pickSoftThreshold(transcriptome_scaled, powerVector = powers,
                         networkType = "signed")
soft_power <- sft$powerEstimate   # 自动选择的软阈值

# 一步法构建网络（适合初学者）
net <- blockwiseModules(
  transcriptome_scaled,
  power = soft_power,
  TOMType = "signed",          # 有符号TOM（保留表达方向）
  minModuleSize = 30,          # 最小模块基因数
  mergeCutHeight = 0.25,       # 模块合并阈值
  numericLabels = FALSE,       # 用颜色命名模块
  verbose = 3
)

# 提取模块特征基因（Module Eigengene）
MEs <- moduleEigengenes(transcriptome_scaled, net$colors)$eigengenes

# 步骤2：模块特征基因 × 菌群丰度相关分析
microbiome_clr <- as.data.frame(clr(microbiome + 0.5))
cor_result <- cor(MEs, microbiome_clr, method = "spearman")  # Spearman相关

# 可视化热图
library(pheatmap)
pheatmap(cor_result,
  main = "基因模块-菌群相关热图",
  fontsize_row = 8,
  fontsize_col = 6
)