mofaplus-multi-omics

MOFA+ 多组学因子分析

概述

MOFA+（多组学因子分析 v2）是一种无监督统计框架，可将多个组学数据集联合分解为一组少量潜在因子。每个因子都捕获一个独立的变异来源（例如细胞周期、疾病表型或技术批次），并与特征权重（载荷）相关联，以揭示哪些基因、峰或蛋白质驱动该因子。Python 包 mofapy2 会生成 HDF5 模型文件，可在 Python 和 R 中用于后续分析。MOFA+ 扩展了原版 MOFA，支持多组设置，其中样本属于不同队列或条件。

适用场景

• 整合来自同一组细胞或样本的两个或多个组学层（例如 scRNA-seq + scATAC-seq、RNA + 蛋白质组、甲基化 + RNA）
• 在没有监督标签的情况下，识别跨组学模态共享的变异来源以及模态特有的变异来源
• 在多组分析中比较不同患者组、治疗条件或时间点的潜在因子差异
• 在聚类、轨迹推断或生存建模前降低多组学数据的维度
• 通过稀疏载荷发现驱动各因子的基因组特征（基因、峰、蛋白质）
• 将因子得分与样本元数据（年龄、分期、治疗反应）相关联，为潜在因子添加注释
• 训练后使用 omics-plotting SKILL，绘制适合发表的因子散点图和载荷热图
• 如果需要借助 VAE 架构，在不同模态间进行深度生成式批次校正并显式推断潜在空间，请改用 scVI / MultiVI（scverse）
• 如果主要目标是通过 iNMF 整合不同技术的数据集（例如 snRNA-seq + snATAC-seq），并得到共享因子和数据集特有因子，请改用 LIGER

前置条件

• Python 包：mofapy2>=0.7、anndata>=0.10、numpy>=1.24、pandas>=2.0、scipy>=1.11、matplotlib>=3.7、seaborn>=0.13、muon（可选，用于 MuData 整合）
• 数据要求：每个组学视图对应一个 AnnData 对象（细胞/样本 x 特征）。所有视图必须共享相同的 obs（细胞/样本）索引。支持各组中存在缺失样本。
• 环境：Python 3.9+；训练好的模型以 HDF5 格式保存，也可通过 MOFA2 Bioconductor 包在 R 中分析

pip install mofapy2 anndata muon matplotlib seaborn

预先访谈

编写任何分析代码之前，请先与用户确认以下事项。

decisions:
  - id: D1
    param: views
    kind: required
    source: data
    ask: "Which omics layers go into the model, and do they share the same samples?"
    default: null

  - id: D2
    param: likelihoods
    kind: required
    source: data
    depends_on: [D1]
    ask: "For each layer, is the data continuous and normalized, binary, or raw counts?"
    default: null

  - id: D3
    param: scaleViews
    kind: required
    source: user
    depends_on: [D1]
    ask: "Should layers be rescaled so the one with the largest variance does not dominate every factor?"
    default: "rescaled"

  - id: D4
    param: numFactors
    kind: required
    source: user
    ask: "How many latent factors should the model start from? Unused ones are pruned automatically."
    default: 15

  - id: D5
    param: groups
    kind: optional
    source: data
    ask: "Should samples be split into groups so the model can compare variance across them?"
    default: "one group"

  - id: D6
    param: convergence
    kind: optional
    source: user
    ask: "Stop early for a quick look, or run to a tight tolerance for the final model?"
    default: "medium tolerance, up to 1000 iterations"

  - id: D7
    param: factorDropThreshold
    kind: optional
    source: user
    ask: "Drop factors that explain almost no variance in any layer?"
    default: "drop below 1% variance explained"

  - id: D8
    param: sparsityPriors
    kind: optional
    source: user
    ask: "Should MOFA retain spike-and-slab weights and per-view relevance priors for sparse, interpretable factors?"
    default: "both enabled"

D2 依赖 D1，而这项决策一旦出错，就会悄无声息地破坏模型：将原始计数交给声明为 gaussian 的层，会使模型拟合测序深度，最终得到的因子反映的是深度，而不是生物学。

快速开始

import numpy as np
import pandas as pd
import anndata as ad
from mofapy2.run.entry_point import entry_point

# Simulate two omics views, 200 cells, 500 RNA genes, 300 ATAC peaks
np.random.seed(42)
n_cells, n_rna, n_atac = 200, 500, 300

adata_rna  = ad.AnnData(np.abs(np.random.randn(n_cells, n_rna)),
                         obs=pd.DataFrame(index=[f"cell_{i}" for i in range(n_cells)]))
adata_atac = ad.AnnData(np.abs(np.random.randn(n_cells, n_atac)),
                         obs=adata_rna.obs.copy())

ent = entry_point()
ent.set_data_options(scale_groups=False, scale_views=True)
ent.set_data_matrix([[adata_rna.X, adata_atac.X]],
                    likelihoods=["gaussian", "gaussian"],
                    views_names=["RNA", "ATAC"],
                    groups_names=["all_cells"],
                    samples_names=[list(adata_rna.obs_names)])
ent.set_model_options(factors=10)
ent.set_train_options(iter=500, convergence_mode="fast", seed=42)
ent.build()
ent.run()
ent.save("mofa_model.hdf5")
print("Model saved to mofa_model.hdf5")

工作流

第 1 步：加载并准备多组学数据

每个组学层都由一个 AnnData 对象表示。对齐不同模态的细胞索引，对 RNA 计数进行对数归一化，并视情况对 ATAC/甲基化数据进行二值化或归一化。

import numpy as np
import pandas as pd
import anndata as ad
import scipy.sparse as sp

# --- RNA-seq: 200 cells x 2000 highly variable genes ---
np.random.seed(42)
n_cells = 200
cell_ids = [f"cell_{i:03d}" for i in range(n_cells)]

# Simulate log-normalized counts (in practice: load from h5ad after Scanpy preprocessing)
rna_counts = np.abs(np.random.randn(n_cells, 2000) * 2)
adata_rna = ad.AnnData(
    X=rna_counts,
    obs=pd.DataFrame(
        {"condition": ["A"] * 100 + ["B"] * 100,
         "patient": [f"P{i % 10}" for i in range(n_cells)]},
        index=cell_ids
    ),
    var=pd.DataFrame(index=[f"Gene_{i}" for i in range(2000)])
)

# --- ATAC-seq: same cells x 1000 peaks ---
atac_matrix = (np.random.rand(n_cells, 1000) > 0.8).astype(float)
adata_atac = ad.AnnData(
    X=atac_matrix,
    obs=adata_rna.obs.copy(),
    var=pd.DataFrame(index=[f"Peak_{i}" for i in range(1000)])
)

# Confirm alignment
assert list(adata_rna.obs_names) == list(adata_atac.obs_names), "Cell indices must match"
print(f"RNA: {adata_rna.shape}, ATAC: {adata_atac.shape}")
print(f"Conditions: {adata_rna.obs['condition'].value_counts().to_dict()}")

第 2 步：创建 MOFA+ 模型对象

实例化 entry_point 并注册所有数据视图。视图按列表的列表形式提供：data[groups][views]。为视图和组分配有意义的名称，以便于解读。

from mofapy2.run.entry_point import entry_point

ent = entry_point()

# Configure data options before setting data
ent.set_data_options(
    scale_groups=False,   # Do not rescale variance between groups
    scale_views=True,     # Rescale each view to unit variance (recommended when views differ in scale)
)

# Provide data as list-of-lists: [groups][views]
# Single group → wrap each view in a list
ent.set_data_matrix(
    data=[[adata_rna.X, adata_atac.X]],       # outer list = groups, inner = views
    likelihoods=["gaussian", "bernoulli"],      # gaussian for continuous, bernoulli for binary ATAC
    views_names=["RNA", "ATAC"],
    groups_names=["all_cells"],
    samples_names=[list(adata_rna.obs_names)]  # one list per group
)

print("Data registered. Views: RNA, ATAC | Groups: all_cells")

第 3 步：设置模型和训练选项

配置因子数量和训练超参数。因子数量越多，捕获的细微变异越多，但计算量和过拟合风险也会增加；convergence_mode="medium" 可兼顾速度和准确性。

# Model options
ent.set_model_options(
    factors=15,          # Number of latent factors (start with 15; prune inactive ones automatically)
    spikeslab_weights=True,  # Sparse weight prior (ARD+spike-slab); recommended for feature selection
    ard_factors=True,    # Automatic relevance determination per factor per view
    ard_weights=True,    # ARD per feature weight; enables pruning of irrelevant features
)

# Training options
ent.set_train_options(
    iter=1000,                     # Maximum EM iterations
    convergence_mode="medium",     # "fast" (<1000 iter), "medium" (default), "slow" (>5000 iter)
    startELBO=1,                   # Start computing ELBO from iteration 1
    freqELBO=5,                    # Compute ELBO every 5 iterations
    dropR2=0.01,                   # Drop factors explaining < 1% variance (set to None to disable)
    seed=42,
    verbose=False
)

print("Model and training options set")

第 4 步：构建并训练模型

构建内部数据结构，然后运行变分推断。训练会生成拟合后的模型，其中每个因子的权重和得分都经过优化，以最大化证据下界（ELBO）。

# Build internal model structure
ent.build()

# Run training (EM algorithm with variational Bayes updates)
ent.run()

# Save trained model to HDF5 — required for downstream analysis
output_path = "mofa_model.hdf5"
ent.save(output_path, overwrite=True)
print(f"Model trained and saved to {output_path}")

第 5 步：加载训练好的模型并检查解释方差

加载 HDF5 模型并检查每个视图中各因子解释的方差。解释总方差低于约 1-2% 的因子通常是噪声。

import h5py
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

def load_mofa_r2(model_path):
    """Extract variance explained (R2) per factor per view from MOFA+ HDF5."""
    with h5py.File(model_path, "r") as f:
        r2 = f["variance_explained"]["r2_per_factor"]
        views = [v.decode() for v in f["views"]["views"][:]]
        groups = [g.decode() for g in f["groups"]["groups"][:]]

        # r2_per_factor: shape (n_groups, n_views, n_factors)
        r2_array = np.stack([r2[g][:] for g in groups], axis=0)  # (groups, views, factors)

    # Average across groups; shape → (n_views, n_factors)
    r2_mean = r2_array.mean(axis=0)
    n_factors = r2_mean.shape[1]
    df = pd.DataFrame(r2_mean * 100,
                      index=views,
                      columns=[f"Factor{i+1}" for i in range(n_factors)])
    return df

r2_df = load_mofa_r2("mofa_model.hdf5")
print("Variance explained (%) per factor per view:")
print(r2_df.round(2))

# Heatmap of variance explained
fig, ax = plt.subplots(figsize=(max(8, r2_df.shape[1] * 0.6), 3))
sns.heatmap(r2_df, annot=True, fmt=".1f", cmap="YlOrRd",
            linewidths=0.5, ax=ax, vmin=0)
ax.set_title("MOFA+ Variance Explained (%) per Factor per View")
ax.set_xlabel("Factor")
ax.set_ylabel("Omics View")
plt.tight_layout()
plt.savefig("mofa_variance_explained.png", dpi=200)
print("Saved mofa_variance_explained.png")

步骤 6：提取因子得分并与元数据相关联

因子得分（Z 矩阵）是各样本在因子空间中的坐标。将得分与连续或分类元数据相关联，以便从生物学角度注释各因子。

import scipy.stats as stats

def load_mofa_factors(model_path):
    """Load factor scores (Z) from MOFA+ HDF5. Returns DataFrame (samples x factors)."""
    with h5py.File(model_path, "r") as f:
        groups = [g.decode() for g in f["groups"]["groups"][:]]
        factors_list = []
        for g in groups:
            z = f["expectations"]["Z"][g][:]   # shape: (n_factors, n_samples)
            samples = [s.decode() for s in f["samples"][g][:]]
            n_factors = z.shape[0]
            df = pd.DataFrame(z.T, index=samples,
                              columns=[f"Factor{i+1}" for i in range(n_factors)])
            factors_list.append(df)
    return pd.concat(factors_list, axis=0)

factors_df = load_mofa_factors("mofa_model.hdf5")
print(f"Factor scores: {factors_df.shape} (cells x factors)")

# Merge with metadata
meta = adata_rna.obs[["condition", "patient"]].copy()
factors_meta = factors_df.join(meta)

# Point-biserial correlation: factor score vs binary metadata
factor_cols = [c for c in factors_meta.columns if c.startswith("Factor")]
print("\nFactor–Condition correlation (eta-squared approximation):")
for fc in factor_cols[:5]:
    groups_vals = [factors_meta.loc[factors_meta["condition"] == g, fc].values
                   for g in factors_meta["condition"].unique()]
    stat, pval = stats.f_oneway(*groups_vals)
    print(f"  {fc}: F={stat:.2f}, p={pval:.3f}")

步骤 7：可视化因子——散点图和特征热图

在此处提取每个因子的特征权重（载荷），然后阅读 skills/data-visualization/omics-plotting/SKILL.md 并遵循其中的绘图方法，使图表风格保持一致：factors_meta（来自步骤 6）→ 因子散点图（→ figures/mofa_factor_scatter.png）；weights_rna → 对顶部 ± 载荷绘制表达热图（→ figures/mofa_rna_weights_heatmap.png）。

def load_mofa_weights(model_path, view_name):
    """Load feature weights (W) for a specific view. Returns DataFrame (features x factors)."""
    with h5py.File(model_path, "r") as f:
        w = f["expectations"]["W"][view_name][:]   # (n_factors, n_features)
        features = [ft.decode() for ft in f["features"][view_name][:]]
        n_factors = w.shape[0]
        df = pd.DataFrame(w.T, index=features,
                          columns=[f"Factor{i+1}" for i in range(n_factors)])
    return df

weights_rna = load_mofa_weights("mofa_model.hdf5", "RNA")
print(f"RNA weights: {weights_rna.shape} (features x factors)")

# Plot inputs are ready; use the omics-plotting SKILL (`skills/data-visualization/omics-plotting/SKILL.md`) to render:
#   factors_meta -> "factor scatter"  -> figures/mofa_factor_scatter.png
#   weights_rna  -> loading heatmap   -> figures/mofa_rna_weights_heatmap.png

步骤 8：下游分析——按因子得分对细胞聚类并进行富集分析

使用因子得分作为低维嵌入进行聚类，并提取每个因子的高权重基因以进行通路富集分析。

from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
import warnings

# --- Cluster cells using factor scores ---
factor_cols = [c for c in factors_df.columns if c.startswith("Factor")]
X_factors = StandardScaler().fit_transform(factors_df[factor_cols].values)

kmeans = KMeans(n_clusters=4, random_state=42, n_init=10)
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    cluster_labels = kmeans.fit_predict(X_factors)

factors_df["mofa_cluster"] = cluster_labels.astype(str)
print(f"K-means clusters (k=4):\n{pd.Series(cluster_labels).value_counts().sort_index()}")

# --- Extract top-weighted genes per factor for enrichment input ---
weights_rna = load_mofa_weights("mofa_model.hdf5", "RNA")

enrichment_inputs = {}
for fc in [f"Factor{i+1}" for i in range(min(5, weights_rna.shape[1]))]:
    # Positive weights: top activating genes; negative: top repressing genes
    top_pos = weights_rna[fc].nlargest(100).index.tolist()
    top_neg = weights_rna[fc].nsmallest(100).index.tolist()
    enrichment_inputs[f"{fc}_positive"] = top_pos
    enrichment_inputs[f"{fc}_negative"] = top_neg

# Save gene lists for external enrichment (e.g., gseapy Enrichr)
for name, genes in list(enrichment_inputs.items())[:2]:
    print(f"\n{name} — top 10: {genes[:10]}")

# Example: run ORA with gseapy (install separately)
# import gseapy
# enr = gseapy.enrichr(gene_list=enrichment_inputs["Factor1_positive"],
#                       gene_sets="GO_Biological_Process_2023",
#                       organism="human", outdir=None)
# print(enr.results.head(5)[["Term", "Adjusted P-value", "Genes"]])

factors_df.to_csv("mofa_factor_scores.csv")
print("\nFactor scores with cluster labels saved to mofa_factor_scores.csv")

关键参数

参数	默认值	范围 / 选项	效果
factors	15	5–50	要推断的潜在因子数量；不活跃的因子会由 ARD 剔除
likelihoods	（必需）	"gaussian"、"bernoulli"、"poisson"	每个视图的似然；归一化连续数据使用 gaussian，二元 ATAC 使用 bernoulli，原始计数使用 poisson
scale_views	True	True/False	将每个视图重新缩放至单位方差；视图的尺度或单位不同时建议启用
scale_groups	False	True/False	对组间方差进行重新缩放；如果各组的总方差差异很大，请设为 True
spikeslab_weights	True	True/False	对权重使用 Spike-and-slab 稀疏先验；通过接近零的权重实现特征选择
ard_factors	True	True/False	对每个视图中的每个因子进行自动相关性判定；剔除未在某个视图中使用的因子
iter	1000	200–5000	EM 迭代的最大次数；通常在 200–800 次迭代内收敛
convergence_mode	"medium"	"fast"、"medium"、"slow"	ELBO 收敛容差：fast = 1e-4，medium = 1e-6，slow = 1e-8
dropR2	0.01	None、0.001–0.05	剔除解释方差低于此比例的因子；None 表示保留所有因子
startELBO	1	1–100	开始监控 ELBO 的迭代次数；设为较高值可跳过初始的不稳定阶段

关键概念

潜在因子和载荷

每个潜在因子 Z_k（长度为 n_samples 的向量）都代表一种变异来源。对应的权重矩阵 W_k（每个视图中长度为 n_features 的向量）包含各特征在该因子上的载荷。使用 Spike-and-slab 先验的因子会产生稀疏载荷：大多数权重会被压缩至接近零，仅保留一小组能够显著驱动该因子的特征。正权重表示因子得分越高，特征表达量也越高。

多组与单组

MOFA+ 中的“组”是一组共享相同因子权重矩阵 W、但因子得分分布 Z 相互独立的样本。以下情况适合使用多组分析：

• 样本来自具有批次级差异的不同队列（患者、数据集）
• 你希望比较各因子在组内与组间的解释程度
• 你预计相同的生物学程序在不同条件下的作用方式不同

如果样本来自没有显著批次结构的单次实验，则适合使用单组分析（所有样本归为一组）。

方差解释热图解读

R2 热图（视图 × 因子）是主要诊断输出。因子应呈现以下特征：

• 视图特异性 R2：仅驱动 RNA 的因子捕捉转录程序；同时驱动 RNA 和 ATAC 的因子捕捉与染色质可及性相关的基因表达
• R2 逐渐降低：因子 1 解释的方差最多；应按顺序检查各因子
• 在所有视图中 R2 均低于 1% 的因子通常可以忽略

常用方法

配方：跨条件的多组分析

适用于比较两个或更多患者组或实验条件，并希望识别条件特异性因素与共享因素的情况。

from mofapy2.run.entry_point import entry_point
import numpy as np
import pandas as pd
import anndata as ad

# Simulate two groups: condition A (100 cells) and condition B (100 cells)
np.random.seed(0)
n_per_group, n_genes, n_peaks = 100, 1000, 500
groups = {"condA": {}, "condB": {}}

for g in groups:
    groups[g]["rna"] = np.abs(np.random.randn(n_per_group, n_genes))
    groups[g]["atac"] = (np.random.rand(n_per_group, n_peaks) > 0.75).astype(float)

sample_ids_A = [f"A_cell_{i}" for i in range(n_per_group)]
sample_ids_B = [f"B_cell_{i}" for i in range(n_per_group)]

ent = entry_point()
ent.set_data_options(scale_groups=False, scale_views=True)

# Multi-group: data[groups][views] — two groups, each with RNA and ATAC
ent.set_data_matrix(
    data=[[groups["condA"]["rna"], groups["condA"]["atac"]],
          [groups["condB"]["rna"], groups["condB"]["atac"]]],
    likelihoods=["gaussian", "bernoulli"],
    views_names=["RNA", "ATAC"],
    groups_names=["condA", "condB"],
    samples_names=[sample_ids_A, sample_ids_B]
)

ent.set_model_options(factors=10, spikeslab_weights=True, ard_factors=True)
ent.set_train_options(iter=500, convergence_mode="fast", seed=0, verbose=False)
ent.build()
ent.run()
ent.save("mofa_multigroup.hdf5", overwrite=True)

# Compare factor scores between groups
factors_all = load_mofa_factors("mofa_multigroup.hdf5")
factors_all["group"] = ["condA"] * n_per_group + ["condB"] * n_per_group
print(f"Multi-group model trained. Factor scores: {factors_all.shape}")
print(factors_all.groupby("group")[["Factor1", "Factor2"]].mean().round(3))

配方：按权重最高的基因识别并注释因素

检索每个因素中权重最高的正向和负向基因，并打印摘要表，以便进行生物学注释。

import pandas as pd

def annotate_factors(model_path, view_name="RNA", top_n=20, n_factors=5):
    """
    Summarize top-loading features per factor to assist biological annotation.

    Returns a DataFrame with factor names, top positive genes, top negative genes,
    and the absolute weight range (an activity proxy).
    """
    weights = load_mofa_weights(model_path, view_name)
    factor_cols = [f"Factor{i+1}" for i in range(min(n_factors, weights.shape[1]))]

    rows = []
    for fc in factor_cols:
        w = weights[fc]
        top_pos = w.nlargest(top_n).index.tolist()
        top_neg = w.nsmallest(top_n).index.tolist()
        rows.append({
            "Factor": fc,
            "Max_weight": round(w.max(), 4),
            "Min_weight": round(w.min(), 4),
            "Top_positive": ", ".join(top_pos[:5]),
            "Top_negative": ", ".join(top_neg[:5]),
        })

    summary = pd.DataFrame(rows)
    return summary

summary_df = annotate_factors("mofa_model.hdf5", view_name="RNA", top_n=20, n_factors=5)
print("\nFactor annotation summary (RNA view):")
print(summary_df.to_string(index=False))
summary_df.to_csv("mofa_factor_annotation.csv", index=False)
print("\nSaved mofa_factor_annotation.csv")

预期输出

输出	描述
mofa_model.hdf5	训练好的 MOFA+ 模型——因子得分、权重、ELBO 轨迹、解释的方差
mofa_variance_explained.png	每个视图中各因子的 R2 (%) 热图；用于选择因子的主要诊断图
figures/mofa_factor_scatter.png	Factor1 与 Factor2/3 的散点图，按元数据（条件、患者）着色；由 omics-plotting 提供
figures/mofa_rna_weights_heatmap.png	前 5 个因子中权重最高的 RNA 特征热图；由 omics-plotting 提供
mofa_factor_scores.csv	每个细胞的因子得分表（细胞数 x 因子数），包含聚类标签
mofa_factor_annotation.csv	因子注释表：每个因子中权重最高的正向/负向基因
各因子的基因列表	可作为 gseapy Enrichr 或 GSEA 的输入，用于识别各因子中富集的通路

故障排除

问题	原因	解决方案
set_data_matrix 中出现 KeyError	组数、视图数或样本列表维度不匹配	确保 data、likelihoods、views_names、groups_names、samples_names 的长度一致
所有因子解释的方差都接近零	数据未预处理或不同视图的尺度不匹配	输入前对每个视图进行归一化；设置 scale_views=True；确认输入矩阵的方差非零
模型完成训练，但因子得分为 NaN	极端值或近奇异数据导致收敛失败	检查输入矩阵中是否有 Inf/NaN；减小 iter；先尝试 convergence_mode="fast"
剪枝掉太多因子（只剩 1-2 个）	dropR2 阈值过于激进，或数据变异不足	设置 dropR2=0.001 或 dropR2=None；提高数据多样性或降低噪声
无法读取 HDF5 文件	训练期间崩溃导致文件被截断	重新运行训练；检查磁盘空间；在 ent.save() 中使用 overwrite=True
所有样本的因子得分都相同	组中只有一个样本，或输入视图方差为零	确认每组至少有 2 个不同样本；检查输入矩阵是否全为零
训练非常缓慢（超过 1 小时）	特征空间过大（每个视图超过 10k 个特征）或因子过多	每个视图预筛选至权重最高的 HVG（2000-5000 个）；将因子数降至 10-15；启用 verbose=False
ELBO 未收敛（持续振荡）	学习率不稳定或数据尺度不佳	增大 startELBO；分别对每个视图进行标准化；使用 convergence_mode="slow"
一个视图的权重都接近零	对连续数据使用了 Bernoulli 似然，或反之	确认 likelihoods 列表与视图数据类型匹配；对归一化后的 RNA 使用 "gaussian"
ModuleNotFoundError: mofapy2	未安装软件包	pip install mofapy2

参考资料

• MOFA2 GitHub 仓库 — mofapy2（Python）和 MOFA2 R 软件包的源代码
• MOFA2 文档和教程 — 官方文档，包含 Python 和 R 的教程
• Argelaguet et al. (2020) — “MOFA+：用于全面整合多模态单细胞数据的统计框架”，Genome Biology
• Argelaguet et al. (2018) — MOFA 原始论文，Molecular Systems Biology
• mofapy2 PyPI 软件包 — 安装说明和版本历史