402_单细胞多组学整合MOFA+¶

# 安装：pip install mofapy2
# 或 R: BiocManager::install("MOFA2")

# ========== Python接口示例 ==========
from mofapy2.run.entry_point import entry_point
import numpy as np
import pandas as pd

# ========== 1. 准备输入数据 ==========
np.random.seed(42)
n_samples = 50   # 样本数
n_genes = 500    # 基因数
n_peaks = 300    # ATAC peak数
n_proteins = 100 # 蛋白数

# 模拟三种组学数据（样本 x 特征）
rna_data = pd.DataFrame(
    np.random.negative_binomial(10, 0.5, size=(n_samples, n_genes)),
    columns=[f"Gene_{i}" for i in range(n_genes)]
)

atac_data = pd.DataFrame(
    np.random.binomial(1, 0.3, size=(n_samples, n_peaks)),  # 0/1二值化数据
    columns=[f"Peak_{i}" for i in range(n_peaks)]
)

protein_data = pd.DataFrame(
    np.random.normal(0, 1, size=(n_samples, n_proteins)),
    columns=[f"Protein_{i}" for i in range(n_proteins)]
)

print(f"RNA数据: {rna_data.shape}")      # (50, 500)
print(f"ATAC数据: {atac_data.shape}")    # (50, 300)
print(f"蛋白数据: {protein_data.shape}") # (50, 100)

# ========== 2. 初始化MOFA+模型 ==========
ent = entry_point()  # 创建MOFA+入口

# 设置数据（字典格式：{组学名: {分组名: 数据矩阵}}）
data_dict = {
    "RNA": {"group1": rna_data.values.T},     # 注意：特征 x 样本
    "ATAC": {"group1": atac_data.values.T},
    "Protein": {"group1": protein_data.values.T}
}
ent.set_data_options(
    scale_groups=False,  # 不对组间缩放
    scale_views=True     # 对不同组学缩放（推荐，因量级不同）
)
ent.set_data_matrix(data_dict)

# ========== 3. 设置模型参数 ==========
ent.set_model_options(
    factors=10,           # 潜在因子数量（从10开始，根据方差解释量调整）
    likelihoods=["gaussian", "bernoulli", "gaussian"]  # 对应RNA/ATAC/蛋白的分布
)

ent.set_train_options(
    iter=1000,      # 最大迭代次数
    convergence_mode="fast",  # fast/medium/slow
    seed=42         # 随机种子（保证可重复）
)

# ========== 4. 训练模型 ==========
ent.build()
ent.run()

# 保存模型到HDF5文件
ent.save("mofa_model.hdf5")
print("模型已保存！")

# ========== 5. 下游分析（R语言更友好） ==========
# R代码示例：
r_code = """
library(MOFA2)  # 加载MOFA2包

# 加载保存的模型
model <- load_model("mofa_model.hdf5")

# 查看每个因子解释的方差比例
plot_variance_explained(model, max_r2=15)  # 热图展示

# 查看因子值（每个样本的因子得分）
factors <- get_factors(model, as.data.frame=TRUE)

# 查看因子权重（每个特征对因子的贡献）
weights <- get_weights(model, as.data.frame=TRUE)

# 找出Factor1权重最大的基因（=对Factor1贡献最大的基因）
top_genes <- weights %>%
    filter(view == "RNA", factor == "Factor1") %>%
    arrange(desc(abs(value))) %>%
    head(20)

# 关联因子与临床信息（如病人分组）
plot_factor(model, factor=1, color_by="condition")
"""
print("R分析代码（复制到R中运行）：")
print(r_code)