多组学整合DIABLO¶

单组学的局限：
  - 转录组只看mRNA，不等于蛋白和代谢物水平
  - 蛋白组遗漏翻译后修饰的调控
  - 代谢组不知道上游基因调控

多组学整合的优势：
  1. 更全面的生物学图景
  2. 跨层次验证（基因→蛋白→代谢物一致变化更可靠）
  3. 发现跨组学的调控关系
  4. 更准确的生物标志物（多组学签名）

三种整合策略：
  ① 早期整合（Early Integration）：拼接数据矩阵后分析
     简单但忽略了组学间差异
  ② 晚期整合（Late Integration）：各组学分别分析后合并结果
     保留组学特性但不能发现跨组学关系
  ③ 中间整合（Intermediate Integration）：DIABLO属于此类
     同时分析多组学数据，发现跨组学相关特征

mixOmics包含的主要方法：

| 方法 | 功能 | 有/无监督 | 数据块数 |
|------|------|-----------|----------|
| PCA | 主成分分析 | 无监督 | 1 |
| PLS | 偏最小二乘 | 有监督 | 2 |
| sPLS | 稀疏PLS | 有监督+变量选择 | 2 |
| PLS-DA | PLS判别分析 | 有监督（分类） | 1+标签 |
| sPLS-DA | 稀疏PLS-DA | 有监督+变量选择 | 1+标签 |
| DIABLO | 多模块sPLS-DA | 有监督+变量选择 | N+标签 |
| MINT | 批次效应校正PLS-DA | 有监督 | 1+标签+批次 |

# 安装mixOmics
BiocManager::install("mixOmics")  # 安装
library(mixOmics)                  # 加载

# DIABLO需要的输入数据：
# 1. 多个组学数据矩阵（行=样本，列=特征）
# 2. 样本分组标签（Y）
# 3. 设计矩阵（控制组学间关联）

# 示例：3个组学数据块
# X_mrna: 100个样本 × 2000个基因
# X_mirna: 100个样本 × 300个miRNA
# X_protein: 100个样本 × 500个蛋白

# 准备数据列表
X <- list(
  mRNA = X_mrna,         # 转录组矩阵
  miRNA = X_mirna,       # miRNA矩阵
  protein = X_protein    # 蛋白组矩阵
)

# 分组标签
Y <- factor(sample_groups)  # 如c("Cancer", "Normal", "Cancer", ...)
cat("样本分组:", table(Y), "\n")

# 检查数据一致性
# 所有矩阵的行名（样本名）必须一致！
lapply(X, dim)  # 查看各矩阵维度
sapply(X, function(x) identical(rownames(x), rownames(X[[1]])))  # 检查样本名一致

# 设计矩阵控制不同组学数据之间的关联强度
# 值为0-1之间：0=完全独立，1=完全关联

# 方法1：全关联设计（推荐起始）
design <- matrix(0.1, nrow = length(X), ncol = length(X),
                 dimnames = list(names(X), names(X)))  # 初始化
diag(design) <- 0  # 对角线为0（自己和自己不需要关联）
design
#         mRNA miRNA protein
# mRNA     0.0   0.1     0.1
# miRNA    0.1   0.0     0.1
# protein  0.1   0.1     0.0

# 方法2：基于先验知识的设计
# 如果知道miRNA调控mRNA但与蛋白关系较远
design_custom <- matrix(c(
  0, 0.8, 0.3,    # mRNA和miRNA强关联，和蛋白弱关联
  0.8, 0, 0.1,    # miRNA和mRNA强关联，和蛋白弱关联
  0.3, 0.1, 0     # 蛋白和mRNA中等关联，和miRNA弱关联
), nrow = 3, dimnames = list(names(X), names(X)))

# 设计矩阵对结果影响：
# 值越大 → 强调组学间的相关性 → 选出跨组学高相关的变量
# 值越小 → 强调判别能力 → 选出区分分组最好的变量
# 通常用0.1作为起始值，然后根据结果调整

# ===== 步骤1：确定每个组学的最佳成分数 =====
# 先用PLS-DA对每个组学单独评估
perf_plsda <- list()
for (i in names(X)) {
  plsda_model <- plsda(X[[i]], Y, ncomp = 10)  # 拟合PLS-DA
  perf_plsda[[i]] <- perf(plsda_model,          # 性能评估
                           validation = "Mval",   # M折交叉验证
                           folds = 5,             # 5折
                           nrepeat = 10)           # 重复10次
}
# 查看每个组学的最佳成分数
# 通常2-5个成分就够了

# ===== 步骤2：确定每个组学每个成分选多少变量 =====
# 使用tune.block.splsda函数
test_keepX <- list(
  mRNA = c(10, 20, 50, 100),     # mRNA候选变量数
  miRNA = c(5, 10, 15, 20),      # miRNA候选变量数
  protein = c(10, 20, 30, 50)    # protein候选变量数
)

tune_result <- tune.block.splsda(
  X = X,                          # 多组学数据
  Y = Y,                          # 分组标签
  ncomp = 2,                      # 成分数
  test.keepX = test_keepX,        # 候选变量数
  design = design,                # 设计矩阵
  validation = "Mval",            # M折交叉验证
  folds = 5,                      # 5折
  nrepeat = 10,                   # 重复10次
  dist = "centroids.dist"         # 距离度量
)

# 查看最佳变量数
tune_result$choice.keepX  # 每个组学每个成分的最佳变量数
# 例如：mRNA=50, miRNA=10, protein=20

# 查看分类错误率
plot(tune_result)  # 画错误率曲线

# 使用最优参数运行DIABLO
keepX <- tune_result$choice.keepX  # 最优变量数
ncomp <- tune_result$choice.ncomp  # 最优成分数

# 运行DIABLO（block.splsda）
diablo_model <- block.splsda(
  X = X,               # 多组学数据列表
  Y = Y,               # 分组标签
  ncomp = ncomp,        # 成分数
  keepX = keepX,        # 每个组学的变量数
  design = design       # 设计矩阵
)

# 查看模型摘要
diablo_model

# 交叉验证评估
perf_diablo <- perf(
  diablo_model,
  validation = "Mval",    # M折交叉验证
  folds = 5,              # 5折
  nrepeat = 10,           # 重复10次
  dist = "centroids.dist" # 距离度量
)

# 查看分类错误率
perf_diablo$error.rate     # 每个组学的错误率
perf_diablo$WeightedVote   # 加权投票的错误率

# 画错误率图
plot(perf_diablo)

# AUC曲线
auc_result <- auroc(diablo_model, roc.block = "miRNA")  # 指定某个组学的AUC

# 1. 个体样本投影图（Score Plot）
plotIndiv(diablo_model,
          ind.names = FALSE,         # 不显示样本名
          legend = TRUE,             # 显示图例
          title = "DIABLO Score Plot",
          comp = c(1, 2))            # 画成分1和2

# 2. 按组学分别画样本图
plotIndiv(diablo_model,
          blocks = "mRNA",           # 只画mRNA组学
          comp = c(1, 2),
          ind.names = FALSE,
          legend = TRUE,
          title = "mRNA Block")

# 1. 变量投影图（Loading Plot）
# 查看每个组学中被选中的变量及其贡献
plotLoadings(diablo_model,
             comp = 1,                # 成分1
             contrib = "max",         # 显示贡献最大的分组
             method = "median")       # 中位数

# 2. 所有组学的变量相关图
plotVar(diablo_model,
        var.names = FALSE,            # 不显示变量名（太多时）
        style = "graphics",
        legend = TRUE,
        pch = c(16, 17, 15),          # 不同组学用不同符号
        cex = c(2, 2, 2))

# 1. Circos图——展示跨组学变量相关性
circosPlot(diablo_model,
           cutoff = 0.7,             # 相关性阈值
           line = TRUE,              # 显示连接线
           color.Y = c("red", "blue"),  # 分组颜色
           size.labels = 1.5)        # 标签大小

# Circos图解读：
# 外圈：不同颜色代表不同组学
# 线条：变量间的相关性（>cutoff才显示）
# 颜色：正相关（红色）/负相关（蓝色）

# 2. 相关网络图
network(diablo_model,
        blocks = c(1, 2, 3),         # 显示所有3个组学
        color.node = c("cyan", "orange", "green"),  # 节点颜色
        cutoff = 0.6,                # 相关性阈值
        save = "network_plot",       # 保存图片
        shape.node = c("circle", "diamond", "square"))  # 节点形状

# 3. 相关热图（CIM图）
cimDiablo(diablo_model,
          color.Y = c("red", "blue"),  # 分组颜色
          margins = c(10, 5),           # 边距
          legend.position = "topright") # 图例位置

# 提取每个组学选中的变量（多组学签名）
selected_vars <- selectVar(diablo_model, comp = 1)  # 成分1

# mRNA签名
mrna_signature <- selected_vars$mRNA$name     # 基因名
mrna_loading <- selected_vars$mRNA$value      # 载荷值
cat("mRNA签名基因数:", length(mrna_signature), "\n")

# miRNA签名
mirna_signature <- selected_vars$miRNA$name
cat("miRNA签名数:", length(mirna_signature), "\n")

# 蛋白签名
protein_signature <- selected_vars$protein$name
cat("蛋白签名数:", length(protein_signature), "\n")

# 导出签名
signature_df <- data.frame(
  omics = c(rep("mRNA", length(mrna_signature)),
            rep("miRNA", length(mirna_signature)),
            rep("protein", length(protein_signature))),
  feature = c(mrna_signature, mirna_signature, protein_signature)
)
write.csv(signature_df, "DIABLO_signature.csv", row.names = FALSE)

# 用训练好的DIABLO模型预测新样本
# 新样本也需要有相同的多组学数据

new_X <- list(
  mRNA = new_mrna_data,      # 新样本mRNA数据
  miRNA = new_mirna_data,    # 新样本miRNA数据
  protein = new_protein_data # 新样本蛋白数据
)

# 预测
predictions <- predict(diablo_model, newdata = new_X)

# 获取预测标签
predicted_class <- predictions$WeightedVote$centroids.dist[, 2]  # 成分2的加权预测
cat("预测结果:", table(predicted_class), "\n")

# 加权投票：综合所有组学的预测结果
# 比单组学预测更准确

# DIABLO分析流程
准备数据列表(X) + 标签(Y) + 设计矩阵(design)
→ tune.block.splsda()（调参）
→ block.splsda()（建模）
→ perf()（评估）
→ 可视化（plotIndiv/plotLoadings/circosPlot/cimDiablo）
→ selectVar()（提取签名）

# 核心函数
block.splsda(): DIABLO建模
tune.block.splsda(): 自动调参
perf(): 交叉验证评估
plotIndiv(): 样本散点图
plotLoadings(): 变量贡献图
circosPlot(): 跨组学关联圆形图
cimDiablo(): 聚类热图
network(): 相关网络图
selectVar(): 提取选中变量
predict(): 预测新样本

# 参数选择指南
设计矩阵值: 0.1起始，高值→强调跨组学关联
成分数ncomp: 通常2-5，看交叉验证误差
keepX: tune函数自动确定
交叉验证: Mval（留一法的改进），5-10折，重复10次

# 与其他多组学方法对比
DIABLO: 有监督+变量选择+多组学（mixOmics）
MOFA+: 无监督+多组学潜因子分析（Python/R）
iCluster: 有监督+多组学聚类（R）
SNF: 相似性网络融合（R）