mofaplus-multi-omics について
MOFA+ は、複数のオミクスデータを統合して解析する無教師あり統計フレームワークです。潜在因子と特徴量の重みで共有・ビュー固有の変動を表し、複数グループの比較にも対応します。
主な機能
- 01複数オミクスビューの同時分解
- 02共有・ビュー固有の変動の検出
- 03複数グループや条件の比較
- 04スパースな荷重による特徴量選択
- 05Python と R で解析できる HDF5 モデル
活用シーン
- 01RNA、ATAC、プロテオミクスの統合
- 02患者群や処理条件間の因子比較
- 03クラスタリング前のオミクス次元削減
- 04因子に関連する遺伝子やタンパク質の特定
インストール
コマンドをコピーしてターミナルで実行します。
npx skills add https:/ / github.com/ jaechang-hits/ sciagent-skills --skill mofaplus-multi-omicsインストール後、Claude Code で/mofaplus-multi-omicsを入力して読み込まれたか確認します。
詳細
MOFA+ マルチオミクス因子分析
概要
MOFA+(Multi-Omics Factor Analysis v2)は、複数のオミクスデータセットを少数の潜在因子に同時に分解する、教師なし統計フレームワークです。各因子は、独立した変動要因(例:細胞周期、疾患表現型、技術的バッチ)を捉え、特徴量の重み(ローディング)と関連付けられています。これにより、どの遺伝子、ピーク、タンパク質がその因子を駆動しているかが明らかになります。Python パッケージの mofapy2 は、Python と R の両方で下流解析に対応した HDF5 モデルファイルを生成します。MOFA+ は元の MOFA を拡張し、サンプルが異なるコホートや条件に属するマルチグループ設定をサポートしています。
使用する場面
- 同じ細胞またはサンプル群から得た、2つ以上のオミクスレイヤーを統合する場合(例: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 オブジェクトを1つ用意します(細胞/サンプル x 特徴量)。すべてのビューで同じ obs(細胞/サンプル)インデックスを共有する必要があります。グループごとにサンプルが欠けていても対応できます。
- 環境:Python 3.9 以降。学習済みモデルは HDF5 として保存され、R の
MOFA2Bioconductor パッケージで解析できます
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 に依存しており、誤った判断をするとモデルに気付かれないまま影響します。ガウス分布として宣言されたレイヤーに生カウントを渡すと、モデルは生物学的な変動ではなくシーケンス深度に適合します。その結果得られる因子は、深度を表すものであり、生物学的な変動ではありません。
クイックスタート
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 | 重みにスパイク・アンド・スラブ疎性事前分布を適用します。ほぼゼロの重みによる特徴量選択が可能になります |
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 のベクトル)には、各特徴量がその因子に対して持つローディングが含まれます。スパイク・アンド・スラブ事前分布を持つ因子では、ローディングが疎になります。大半の重みはほぼゼロまで縮小され、因子を実質的に駆動する少数の特徴量だけが残ります。重みが正である場合、因子スコアが高いほど特徴量の発現も高いことを意味します。
マルチグループとシングルグループ
MOFA+ における「グループ」とは、同じ因子重み行列 W を共有しながら、独立した因子スコア分布 Z を持つサンプルの集合です。以下の場合は、マルチグループ解析を使用します。
- サンプルが、バッチレベルの違いがある別個のコホート(患者、データセット)に由来する場合
- 各因子がグループ内とグループ間で説明する分散の量を比較したい場合
- 同じ生物学的プログラムが条件によって異なる形で働くと想定される場合
シングルグループ解析(すべてのサンプルを 1 つのグループにまとめる)は、大きなバッチ構造のない単一の実験からサンプルを得ている場合に適しています。
説明分散ヒートマップの解釈
R2 ヒートマップ(ビュー × 因子)は、主要な診断出力です。因子には、以下の特徴が見られるはずです。
- ビュー固有の R2:RNA のみを駆動する因子は転写プログラムを捉えます。RNA と ATAC の両方を駆動する因子は、クロマチンアクセシビリティと連動した遺伝子発現を捉えます
- R2 の低下:因子 1 が最も多くの分散を説明します。因子は順番に調べてください
- すべてのビューで R2 が 1% 未満の因子は、通常無視できます
よく使われるレシピ
レシピ:条件間の複数グループ解析
2つ以上の患者グループまたは実験条件を比較し、条件固有の要因と共通の要因を特定したい場合に使用します。
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パッケージ — インストール方法とバージョン履歴