751. 多组学单细胞因子分析MOFA+¶

假设存在一些"潜在因子"同时驱动多种组学的变异

Factor 1 可能代表"细胞分化状态"
  → 同时影响分化基因的表达(RNA)
  → 同时影响分化增强子的开放(ATAC)
  → 同时影响分化标记蛋白的表达(蛋白)

Factor 2 可能代表"细胞周期"
  → 同时影响周期基因表达(RNA)
  → 同时影响复制起点的开放(ATAC)

# Python安装
pip install mofapy2  # 安装MOFA+

# R安装
# BiocManager::install("MOFA2")

# ===== MOFA+完整分析流程 =====
from mofapy2.run.entry_point import entry_point  # 导入入口
import pandas as pd
import numpy as np
import scanpy as sc
import matplotlib.pyplot as plt

# ===== 第一步：准备多组学数据 =====
# 假设我们有RNA和ATAC两种数据

# RNA数据
adata_rna = sc.read_h5ad("scRNA.h5ad")
sc.pp.normalize_total(adata_rna, target_sum=1e4)
sc.pp.log1p(adata_rna)
sc.pp.highly_variable_genes(adata_rna, n_top_genes=3000)
rna_matrix = adata_rna[:, adata_rna.var.highly_variable].X.toarray()

# ATAC数据（peaks × cells）
adata_atac = sc.read_h5ad("scATAC.h5ad")
sc.pp.normalize_total(adata_atac, target_sum=1e4)
sc.pp.log1p(adata_atac)
sc.pp.highly_variable_genes(adata_atac, n_top_genes=5000)
atac_matrix = adata_atac[:, adata_atac.var.highly_variable].X.toarray()

# 确保细胞匹配
common_cells = adata_rna.obs_names.intersection(adata_atac.obs_names)
rna_matched = rna_matrix[adata_rna.obs_names.isin(common_cells)]
atac_matched = atac_matrix[adata_atac.obs_names.isin(common_cells)]

# ===== 第二步：创建MOFA+对象 =====
ent = entry_point()  # 创建入口点

# 设置数据
ent.set_data_options(
    scale_groups=False,  # 不跨组缩放
    scale_views=False    # 不跨视图缩放
)

# 添加数据
# MOFA+需要的格式：list of matrices
# 每个视图(view)一个矩阵，每个组(group)一个
data = [[rna_matched],   # 视图1: RNA
        [atac_matched]]   # 视图2: ATAC

ent.set_data_matrix(
    data,
    views_names=["RNA", "ATAC"],   # 视图名称
    groups_names=["all_samples"],   # 组名称
    samples_names=[common_cells.tolist()],  # 样本名
    features_names=[
        adata_rna.var_names[adata_rna.var.highly_variable].tolist(),  # RNA特征名
        adata_atac.var_names[adata_atac.var.highly_variable].tolist()  # ATAC特征名
    ]
)

# ===== 第三步：设置模型参数 =====
ent.set_model_options(
    factors=15,          # 因子数量
    spikeslab_weights=True,  # 使用spike-and-slab先验（促进稀疏性）
    ard_weights=True     # 自动相关性检测（自动去除不重要的因子）
)

ent.set_training_options(
    iter=1000,           # 最大迭代次数
    convergence_mode="fast",  # 收敛模式
    seed=42,             # 随机种子
    gpu_mode=False       # 是否使用GPU
)

# ===== 第四步：训练模型 =====
ent.build()  # 构建模型
ent.run()    # 训练

# 保存模型
ent.save("mofa_model.hdf5")  # 保存为HDF5文件

# ===== R语言分析MOFA+结果（可视化更方便）=====
library(MOFA2)
library(ggplot2)

# 加载训练好的模型
model <- load_model("mofa_model.hdf5")

# 查看模型基本信息
print(model)
# Number of views: 2 (RNA, ATAC)
# Number of groups: 1
# Number of factors: 15
# Number of features: RNA=3000, ATAC=5000

# ===== 1. 方差分解图 =====
# 每个因子在每种组学中解释了多少变异
plot_variance_explained(model, max_r2=15)
# 关键信息：
# Factor 1 在RNA中解释10%变异，在ATAC中解释8% → 共享因子
# Factor 3 只在RNA中解释5% → RNA特异因子
# Factor 5 只在ATAC中解释3% → ATAC特异因子

# ===== 2. 因子相关性 =====
plot_factor_cor(model)
# 因子之间应该不太相关（MOFA+保证正交性）

# ===== 3. 因子可视化（类似UMAP）=====
# 用因子作为坐标
plot_factors(
  model,
  factors = c(1, 2),  # 用前两个因子
  color_by = "cell_type",  # 按细胞类型着色
  dot_size = 1
)

# ===== 4. 权重分析（哪些特征驱动哪个因子）=====
# 查看Factor 1的RNA权重（贡献最大的基因）
plot_weights(
  model,
  view = "RNA",
  factor = 1,
  nfeatures = 20,  # Top 20特征
  scale = TRUE
)

# 查看Factor 1的ATAC权重（贡献最大的peaks）
plot_weights(
  model,
  view = "ATAC",
  factor = 1,
  nfeatures = 20
)

# ===== 5. 特征热图 =====
plot_data_heatmap(
  model,
  view = "RNA",
  factor = 1,
  features = 50,  # Top 50特征
  cluster_rows = TRUE,
  cluster_cols = FALSE,
  show_colnames = FALSE
)

# ===== 6. 因子的GO富集分析 =====
library(clusterProfiler)
library(org.Hs.eg.db)

# 提取Factor 1的top基因
weights <- get_weights(model, views="RNA", factors=1, as.data.frame=TRUE)
top_genes <- weights %>% 
  arrange(desc(abs(value))) %>% 
  head(200) %>% 
  pull(feature)

# GO富集
ego <- enrichGO(
  gene = top_genes,
  OrgDb = org.Hs.eg.db,
  keyType = "SYMBOL",
  ont = "BP",
  pAdjustMethod = "BH",
  pvalueCutoff = 0.05
)
dotplot(ego, showCategory=15)

# ===== MOFA+速查 =====

# Python安装/训练
pip install mofapy2
ent = entry_point()
ent.set_data_matrix(data, views_names=["RNA","ATAC"])
ent.set_model_options(factors=15)
ent.build(); ent.run()
ent.save("model.hdf5")

# R分析结果
library(MOFA2)
model <- load_model("model.hdf5")
plot_variance_explained(model)       # 方差分解
plot_factors(model, factors=c(1,2))  # 因子散点
plot_weights(model, view="RNA", factor=1)  # 权重

# 核心概念：
# View = 组学类型（RNA, ATAC, 蛋白...）
# Group = 样本/条件
# Factor = 潜在变异来源
# Weight = 特征对因子的贡献

# 因子数选择：
# 通常10-20个，看方差解释图
# ARD先验会自动"关掉"不重要的因子