多组学因子分析MOFA¶

# 多组学整合方法分类
# 1. 早期整合（Early integration）：拼接特征矩阵后降维
# 2. 晚期整合（Late integration）：各组学独立分析后合并结果
# 3. 中间整合（Intermediate integration）：联合建模 — MOFA属于此类

# MOFA+模型的数学表达
# Y_m = W_m × Z^T + ε_m
# Y_m: 第m个组学的数据矩阵 (features × samples)
# W_m: 第m个组学的权重矩阵 (features × factors)
# Z:   因子矩阵 (samples × factors)，跨组学共享
# ε_m: 噪声

# === 数据准备示例 ===
library(MOFA2)

# 模拟三组学数据准备
# 1. RNA-seq数据（已为log2 CPM）
rna_data <- read.csv("rna_log2cpm.csv", row.names = 1)
# 取top 5000高变异基因
rna_var <- apply(rna_data, 1, var)
rna_top <- rna_data[names(sort(rna_var, decreasing = TRUE)[1:5000]), ]

# 2. 蛋白组数据（已为log2强度）
prot_data <- read.csv("protein_log2.csv", row.names = 1)
# 去除缺失>50%的蛋白
prot_missing <- rowMeans(is.na(prot_data))
prot_filtered <- prot_data[prot_missing < 0.5, ]

# 3. 甲基化数据（转为M值）
meth_beta <- read.csv("methylation_beta.csv", row.names = 1)
meth_m <- log2(meth_beta / (1 - meth_beta))
# 取top 5000高变异CpG
meth_var <- apply(meth_m, 1, var, na.rm = TRUE)
meth_top <- meth_m[names(sort(meth_var, decreasing = TRUE)[1:5000]), ]

# 确保样本ID一致
common_samples <- Reduce(intersect, list(
  colnames(rna_top),
  colnames(prot_filtered),
  colnames(meth_top)
))

# 中心化（每个特征减去均值）
rna_centered <- t(scale(t(rna_top[, common_samples]), center = TRUE, scale = FALSE))
prot_centered <- t(scale(t(prot_filtered[, common_samples]), center = TRUE, scale = FALSE))
meth_centered <- t(scale(t(meth_top[, common_samples]), center = TRUE, scale = FALSE))

# 创建数据列表
data_list <- list(
  RNA = as.matrix(rna_centered),
  Protein = as.matrix(prot_centered),
  Methylation = as.matrix(meth_centered)
)

# === MOFA+ 模型构建与训练 ===
library(MOFA2)

# 创建MOFA对象
mofa <- create_mofa(data_list)

# 查看数据概况
print(mofa)
plot_data_overview(mofa)

# 设置数据选项
data_opts <- get_default_data_options(mofa)
data_opts$scale_views <- TRUE  # 是否对各组学做方差缩放
# TRUE: 各组学等权重贡献；FALSE: 高方差组学贡献更大

# 设置模型选项
model_opts <- get_default_model_options(mofa)
model_opts$num_factors <- 15  # 初始因子数（会自动删减非活跃因子）
model_opts$likelihoods <- c("gaussian", "gaussian", "gaussian")
# 可选: gaussian（连续）, poisson（count）, bernoulli（二元）

# 设置训练选项
train_opts <- get_default_training_options(mofa)
train_opts$convergence_mode <- "slow"  # slow更精确，fast更快
train_opts$maxiter <- 1000
train_opts$seed <- 42
train_opts$gpu_mode <- FALSE  # 有GPU则设为TRUE

# 训练模型
mofa <- prepare_mofa(
  mofa,
  data_options = data_opts,
  model_options = model_opts,
  training_options = train_opts
)

# 使用Python后端训练（需要mofapy2）
mofa <- run_mofa(mofa, outfile = "mofa_model.hdf5", use_basilisk = TRUE)

# === 模型质量评估 ===

# 方差分解
plot_variance_explained(mofa, max_r2 = 15)
plot_variance_explained(mofa, plot_total = TRUE)  # 各因子总方差

# 查看因子值分布
plot_factor(mofa, factor = 1, color_by = "group_label")

# 因子间相关性（应该接近零）
plot_factor_cor(mofa)

# === 因子权重（Weights）分析 ===

# 各组学中每个因子的top权重特征
plot_weights(mofa, view = "RNA", factor = 1, nfeatures = 20, scale = TRUE)
plot_weights(mofa, view = "Protein", factor = 1, nfeatures = 20, scale = TRUE)

# 权重热图
plot_top_weights(mofa, view = "RNA", factor = 1, nfeatures = 30)

# 提取权重矩阵
weights <- get_weights(mofa, views = "all", factors = "all", as.data.frame = TRUE)

# 提取因子值
factors <- get_factors(mofa, factors = "all", as.data.frame = TRUE)

# === 因子-特征数据热图 ===
plot_data_heatmap(mofa, view = "RNA", factor = 1,
                   features = 50, cluster_rows = TRUE,
                   show_rownames = FALSE, annotation_samples = "group_label")

# === 功能富集分析 ===
library(MOFA2)

# GSEA富集（需要msigdbr或自定义基因集）
library(msigdbr)

# 获取Hallmark基因集
hallmark <- msigdbr(species = "Homo sapiens", category = "H") %>%
  dplyr::select(gs_name, gene_symbol) %>%
  as.data.frame()

# MOFA内置的GSEA
enrichment_results <- run_enrichment(mofa,
  view = "RNA",
  feature.sets = hallmark,
  sign = "all",  # all, positive, negative
  statistical.test = "parametric"
)

# 富集结果可视化
plot_enrichment(enrichment_results, factor = 1, max.pathways = 15)
plot_enrichment_heatmap(enrichment_results)

# === 因子与临床表型关联 ===
# 添加样本元数据
samples_metadata(mofa) <- merge(
  samples_metadata(mofa),
  clinical_data,
  by.x = "sample", by.y = "sample_id"
)

# 因子与连续变量关联
plot_factor(mofa, factor = 1, color_by = "age", add_violin = FALSE,
            dodge = TRUE)

# 因子与分类变量关联
plot_factor(mofa, factor = 1, color_by = "tumor_stage",
            dodge = TRUE, add_violin = TRUE)

# 因子与生存分析
library(survival)
library(survminer)

factors_df <- get_factors(mofa, as.data.frame = TRUE) %>%
  tidyr::spread(factor, value)

surv_data <- merge(factors_df, clinical_data, by.x = "sample", by.y = "sample_id")
surv_data$Factor1_group <- ifelse(surv_data$Factor1 > median(surv_data$Factor1),
                                   "High", "Low")

fit <- survfit(Surv(OS_time, OS_status) ~ Factor1_group, data = surv_data)
ggsurvplot(fit, data = surv_data, pval = TRUE, risk.table = TRUE)

# === MCIA 分析流程 ===
library(omicade4)

# 准备数据列表（每个元素为基因×样本矩阵）
data_mcia <- list(
  RNA = as.data.frame(rna_centered),
  Protein = as.data.frame(prot_centered),
  Methylation = as.data.frame(meth_centered)
)

# 运行MCIA
mcia_result <- mcia(data_mcia, cia.nf = 10)

# 查看各轴解释的方差
plot(mcia_result$mcoa, what = "eigenvalues")

# 样本空间投影（整合视图）
plot(mcia_result, axes = c(1, 2), phenovec = group_labels,
     sample.lab = FALSE)

# 各组学的一致性评估（RV系数）
# RV系数衡量两个数据集之间的全局相似性，值越高越一致
mcia_result$mcoa$RV

# 提取MCIA得分（类似因子值）
mcia_scores <- mcia_result$mcoa$SynVar  # 综合变量得分

# === MOFA+ 多组设计 ===

# 准备多组数据（关键：用groups参数指定分组，而非嵌套列表）
# 假设有两个时间点的多组学数据，样本名分别带有T1/T2前缀
# 方法1：用矩阵列表 + groups向量
data_list_views <- list(
  RNA = rbind(rna_data_T1, rna_data_T2),        # 合并两个时间点
  Protein = rbind(prot_data_T1, prot_data_T2)
)
# groups向量：每个样本属于哪个组
groups <- c(rep("Timepoint1", ncol(rna_data_T1)),
            rep("Timepoint2", ncol(rna_data_T2)))

# 创建多组MOFA对象（通过groups参数传入分组信息）
mofa_multi <- create_mofa(data_list_views, groups = groups)

# 方法2：用长格式data.frame（含group列）
# df_long <- data.frame(sample, feature, view, group, value)
# mofa_multi <- create_mofa(df_long)

# 训练（同上设置）
mofa_multi <- prepare_mofa(mofa_multi,
  data_options = data_opts,
  model_options = model_opts,
  training_options = train_opts)
mofa_multi <- run_mofa(mofa_multi)

# 比较因子在不同组中的方差贡献
plot_variance_explained(mofa_multi, split_by = "group")

# 查看因子值在不同组中的分布
plot_factor(mofa_multi, factor = 1, color_by = "group",
            dodge = TRUE, add_violin = TRUE)

# 注意：compare_factors()用于比较不同MOFA模型间的因子相关性，
# 而非同一模型内的组间比较。组间比较应通过plot_factor()可视化
# 或提取因子值后用统计检验（如Wilcoxon test）进行

# 组特异性因子识别：查看各因子在各组各组学中的方差贡献
plot_variance_explained(mofa_multi, split_by = "group")

# === 安装 ===
# R包
# install.packages("BiocManager")
# BiocManager::install("MOFA2")
# BiocManager::install("omicade4")

# Python后端
pip install mofapy2

# === 常用R函数速查 ===
# create_mofa()                    # 创建MOFA对象
# prepare_mofa()                   # 设置参数
# run_mofa()                       # 训练模型
# plot_data_overview()             # 数据概览
# plot_variance_explained()        # 方差分解
# get_factors() / get_weights()    # 提取结果
# plot_weights() / plot_top_weights()  # 权重可视化
# run_enrichment()                 # GSEA富集分析
# plot_factor() / plot_factors()   # 因子可视化