muon-multiomics-singlecell

Name: muon-multiomics-singlecell
Author: jaechang-hits

github.com/jaechang-hitsデータ・データベース38479 回インストール6 日前に更新

リポジトリ

muon-multiomics-singlecell について

muonは、複数モダリティの単一細胞データを解析する研究者向けのPythonフレームワークです。MuDataコンテナでデータをまとめ、モダリティ別の前処理、WNNによる統合解析、モダリティ間解析を行えます。

主な機能

01MuDataで複数モダリティを格納
02ATACとタンパク質データの前処理
03WNNの統合グラフと埋め込みの構築
04モダリティ間の特徴解析
05scanpy、scvi-tools、MOFA+との連携

活用シーン

01ペアのRNA・ATACデータの解析
02CITE-seqのRNA・タンパク質データの処理
03統合埋め込みと細胞クラスタの作成
04ATACピークと遺伝子発現の関連付け

インストール

コマンドをコピーしてターミナルで実行します。

$npx skills add https://github.com/jaechang-hits/sciagent-skills --skill muon-multiomics-singlecell

インストール後、Claude Code で/muon-multiomics-singlecellを入力して読み込まれたか確認します。

詳細

muon — マルチモーダル単一細胞解析

概要

muon は、AnnData エコシステムを拡張する、マルチモーダル単一細胞データ解析用の Python フレームワークです。中核となるデータ構造 MuData は、複数の AnnData オブジェクト（RNA、ATAC、タンパク質など、モダリティごとに 1 つ）を、共通の観測軸と変数軸を持つ形で格納し、すべてのモダリティにわたる連携した処理を可能にします。muon は、モダリティに応じた前処理ルーチン（ATAC の TF-IDF と LSI、表面タンパク質の CLR 正規化）、統合次元削減のための加重近傍（WNN）グラフ構築、クロスモーダル解析ツールを提供します。scanpy、scvi-tools、MOFA+ と直接連携し、マルチオミクス単一細胞解析の一連のワークフローを実現します。

使用する場面

10x Genomics Multiome データ（同じ核から RNA と ATAC を同時に取得）の解析
CITE-seq 実験（同じ細胞から RNA と表面タンパク質を取得）の処理
WNN を使って 2 つ以上のモダリティのシグナルを統合した、統合 UMAP 埋め込みの構築
ATAC-seq モダリティの前処理（TF-IDF 正規化、LSI 次元削減）
中心化対数比（CLR）正規化による表面タンパク質データの正規化
クロスモーダルな特徴量の対応付け（ATAC ピークと近傍の遺伝子発現の関連付け）
統合コンテナ内の複数のオミクス層に対する MOFA+ 因子分析の適用
一般的な結果テーブルの図には omics-plotting SKILL を使用してください。統合埋め込みには scanpy/muon の sc.pl.* / mu.pl.* を使用します
同時に測定されたオミクスがなく、単一の RNA-seq モダリティを解析する場合は、代わりに scanpy-scrna-seq を使用してください
モダリティ間の確率的な深層生成バッチ補正が必要な場合は、scvi-tools (MultiVI / totalVI) を使用してください

前提条件

Python パッケージ：muon>=0.1.6、scanpy>=1.10、anndata>=0.10、numpy、scipy、pandas、matplotlib、leidenalg
データ要件：10x Multiome h5 または h5mu ファイル、あるいは同じ obs_names を共有する、モダリティごとの AnnData オブジェクト（細胞 x 特徴量）
環境：Python 3.9 以上。50k 細胞を超えるデータセットでは 16 GB 以上の RAM が必要です。MOFA 因子分析用の mofapy2 は任意です

pip install "muon[all]" "scanpy[leiden]" anndata
# Optional: for MOFA+ integration
pip install mofapy2

事前確認

解析コードを書く前に、ユーザーと以下の点を確認してください。

decisions:
  - id: D1
    param: modalities
    kind: required
    source: data
    ask: "Which measurement types are being integrated, and are they from the same cells or matched separately?"
    default: null

  - id: D2
    param: perModalityEmbedding
    kind: required
    source: data
    depends_on: [D1]
    ask: "How should each modality be reduced before joint analysis - PCA for expression, latent semantic indexing for accessibility?"
    default: null

  - id: D3
    param: integrationMethod
    kind: required
    source: user
    depends_on: [D1]
    ask: "Join the modalities by a weighted nearest-neighbour graph, or by a shared factor model?"
    default: "weighted nearest neighbours"

  - id: D4
    param: clusterResolution
    kind: required
    source: user
    ask: "How finely should the joint graph be split into cell groups?"
    default: 1.0

  - id: D5
    param: neighborhoodSize
    kind: optional
    source: user
    ask: "How many neighbours should define each cell's local structure?"
    default: 30

  - id: D6
    param: accessibilityNormalization
    kind: optional
    source: user
    ask: "Which scaling should the accessibility modality use before dimensionality reduction?"
    default: "TF-IDF at 1e4"

  - id: D7
    param: factorCount
    kind: optional_conditional
    source: user
    depends_on: [D3]
    ask: "How many latent factors should the shared factor model infer?"
    default: 10
    skip_if: "weighted nearest-neighbour integration"

D2 は形式的な確認ではありません。統合解析の前にアクセシビリティデータへ発現データ向けの PCA を適用すると、シーケンス深度が第 1 成分を支配するようになり、前処理によって情報量が低くなったモダリティが、情報を持たないという理由で加重グラフによりさらに低く重み付けされます。

クイックスタート

import muon as mu
import scanpy as sc
import numpy as np
import anndata as ad
import pandas as pd

# Simulate a small RNA + ATAC MuData object (200 cells)
np.random.seed(42)
n_cells = 200
rna  = ad.AnnData(np.abs(np.random.negative_binomial(5, 0.3, (n_cells, 2000))).astype(float),
                  obs=pd.DataFrame(index=[f"cell_{i}" for i in range(n_cells)]),
                  var=pd.DataFrame(index=[f"gene_{j}" for j in range(2000)]))
atac = ad.AnnData(np.abs(np.random.negative_binomial(2, 0.5, (n_cells, 5000))).astype(float),
                  obs=rna.obs.copy(),
                  var=pd.DataFrame(index=[f"peak_{j}" for j in range(5000)]))

mdata = mu.MuData({"rna": rna, "atac": atac})
print(mdata)
# MuData object with n_obs × n_vars = 200 × 7000
#   2 modalities
#     rna: 200 × 2000
#     atac: 200 × 5000

# Minimal preprocessing
sc.pp.normalize_total(mdata["rna"], target_sum=1e4)
sc.pp.log1p(mdata["rna"])
sc.pp.highly_variable_genes(mdata["rna"], n_top_genes=500)
sc.pp.pca(mdata["rna"])

mu.atac.pp.tfidf(mdata["atac"])
mu.atac.tl.lsi(mdata["atac"])

# WNN joint embedding
mu.pp.neighbors(mdata, key_added="wnn",
                use_rep={"rna": "X_pca", "atac": "X_lsi"})
sc.tl.umap(mdata, neighbors_key="wnn")
sc.tl.leiden(mdata, neighbors_key="wnn", key_added="leiden_wnn")
print(f"Clusters: {mdata.obs['leiden_wnn'].nunique()}")

コア API

モジュール 1：MuData の作成と入出力

MuData は中核となるコンテナです。複数の AnnData モダリティを格納する辞書と、共有の obs および var スロットで構成されます。mdata.obs の DataFrame には、モダリティごとの観測メタデータが結合されます。曖昧な列にはモダリティのプレフィックスが付与されます。

import muon as mu
import anndata as ad
import numpy as np
import pandas as pd

# --- Build from per-modality AnnData objects ---
n_cells = 300
rna  = ad.AnnData(np.abs(np.random.negative_binomial(5, 0.3, (n_cells, 3000))).astype(float),
                  obs=pd.DataFrame(index=[f"cell_{i}" for i in range(n_cells)]),
                  var=pd.DataFrame(index=[f"gene_{j}" for j in range(3000)]))
prot = ad.AnnData(np.abs(np.random.randn(n_cells, 30)) + 3,
                  obs=rna.obs.copy(),
                  var=pd.DataFrame(index=[f"protein_{j}" for j in range(30)]))

mdata = mu.MuData({"rna": rna, "protein": prot})
print(mdata)
# Access individual modalities
print(mdata.mod["rna"])        # AnnData: 300 × 3000
print(mdata["rna"])            # same shorthand
print(mdata.obs.head())        # shared observation metadata
print(mdata.obs_names[:5])     # cell barcodes

# Save and load
mdata.write("multiome_data.h5mu")
mdata2 = mu.read("multiome_data.h5mu")
print(f"Loaded: {mdata2.n_obs} cells, {mdata2.n_mod} modalities")

# --- Load from 10x Multiome h5 file ---
# mdata = mu.read_10x_h5("filtered_feature_bc_matrix.h5")
# Produces MuData with mdata["rna"] and mdata["atac"] modalities

# --- Subsetting by cells or features ---
# Select high-quality cells (e.g., after QC)
mask = np.ones(mdata.n_obs, dtype=bool)  # replace with actual QC mask
mdata_filtered = mdata[mask].copy()
print(f"After filtering: {mdata_filtered.n_obs} cells")

# Propagate obs mask to each modality
mu.pp.intersect_obs(mdata)   # ensures obs consistency across modalities

モジュール 2：RNA モダリティの前処理

scanpy 関数を使い、mdata["rna"] に標準的な scRNA-seq 前処理を適用します。RNA モダリティは単独の scanpy ワークフローと同様に前処理されますが、MuData コンテナのスライスに対して処理を行います。

import scanpy as sc
import numpy as np

# Assume mdata["rna"] has raw integer counts
rna = mdata["rna"]

# QC metrics
sc.pp.calculate_qc_metrics(rna, percent_top=None, log1p=False, inplace=True)
rna.obs["pct_counts_mt"] = rna[:, rna.var_names.str.startswith("MT-")].X.sum(axis=1).A1 / rna.obs["total_counts"] * 100

# Filter cells and genes
min_genes, max_genes, max_mt = 200, 5000, 20
sc.pp.filter_cells(rna, min_genes=min_genes)
sc.pp.filter_genes(rna, min_cells=5)
rna = rna[(rna.obs["n_genes_by_counts"] < max_genes) &
          (rna.obs["pct_counts_mt"] < max_mt)].copy()
print(f"RNA after QC: {rna.n_obs} cells × {rna.n_vars} genes")

# Normalize and log-transform
sc.pp.normalize_total(rna, target_sum=1e4)
sc.pp.log1p(rna)

# Highly variable genes
sc.pp.highly_variable_genes(rna, n_top_genes=3000, flavor="seurat_v3",
                             batch_key=None)
print(f"HVGs: {rna.var['highly_variable'].sum()}")

# PCA on HVGs
sc.pp.pca(rna, n_comps=50, use_highly_variable=True)
print(f"PCA embedding shape: {rna.obsm['X_pca'].shape}")
# Update slice in MuData
mdata.mod["rna"] = rna

モジュール 3：ATAC モダリティの前処理

ATAC-seq モダリティには、異なる正規化戦略が必要です。TF-IDF（単語頻度–逆文書頻度）は、細胞とピークにわたるピークアクセシビリティを正規化します。LSI（潜在的意味インデックス。TF-IDF の後に適用する切り詰め SVD と同等）は、低次元埋め込みを生成します。通常、第 1 LSI 成分は生物学的特徴ではなくシーケンス深度を捉えるため、除外します。

# Assume mdata["atac"] contains raw binary or integer peak accessibility counts
atac = mdata["atac"]

# Basic QC: filter low-coverage cells and low-frequency peaks
sc.pp.calculate_qc_metrics(atac, percent_top=None, log1p=False, inplace=True)
sc.pp.filter_cells(atac, min_genes=200)     # min peaks detected
sc.pp.filter_genes(atac, min_cells=10)      # min cells a peak appears in
print(f"ATAC after QC: {atac.n_obs} cells × {atac.n_vars} peaks")

# TF-IDF normalization (log(TF) * log(IDF) scaling)
mu.atac.pp.tfidf(atac, scale_factor=1e4)
print("TF-IDF normalization complete")
print(f"ATAC data range: [{atac.X.min():.2f}, {atac.X.max():.2f}]")

# LSI dimensionality reduction (truncated SVD on TF-IDF matrix)
mu.atac.tl.lsi(atac, n_comps=50, use_highly_variable=False)
# LSI component 1 correlates with sequencing depth — exclude it
# mu.atac.tl.lsi sets X_lsi starting from component 2 by default
print(f"LSI embedding shape: {atac.obsm['X_lsi'].shape}")

# Update modality in MuData
mdata.mod["atac"] = atac

モジュール 4：WNN グラフと統合埋め込み

加重近傍（WNN）は、細胞ごと、モダリティごとの重みを学習することで、複数のモダリティの埋め込みを統合します。RNA シグナルの品質が高い細胞では RNA の重みが高くなり、ATAC シグナルがよりクリーンな細胞では ATAC の重みが高くなります。生成された WNN グラフは、UMAP レイアウトと Leiden クラスタリングに使用されます。

# Compute per-modality neighbor graphs first (optional but enables modality-specific UMAPs)
mu.pp.neighbors(mdata["rna"],  use_rep="X_pca",  n_neighbors=30, key_added="neighbors")
mu.pp.neighbors(mdata["atac"], use_rep="X_lsi",  n_neighbors=30, key_added="neighbors")

# WNN joint neighbor graph across modalities
mu.pp.neighbors(
    mdata,
    key_added="wnn",
    use_rep={"rna": "X_pca", "atac": "X_lsi"},
    n_neighbors=30,
    random_state=42,
)
print("WNN graph built. Keys:", list(mdata.obsp.keys()))
# Expected: ['wnn_connectivities', 'wnn_distances']

# UMAP from WNN graph
sc.tl.umap(mdata, neighbors_key="wnn", random_state=42)
print(f"UMAP embedding shape: {mdata.obsm['X_umap'].shape}")

# Leiden clustering from WNN graph
sc.tl.leiden(mdata, neighbors_key="wnn", resolution=0.5, key_added="leiden_wnn")
n_clusters = mdata.obs["leiden_wnn"].nunique()
print(f"Leiden WNN clustering: {n_clusters} clusters at resolution 0.5")

Module 5: 可視化

muon は、モダリティを考慮した関数で scanpy のプロットインターフェースを拡張します。mu.pl.embedding() は、任意のモダリティの特徴量で統合 UMAP 埋め込みを色分けします。また、モダリティの接頭辞を付けた特徴量名（例："rna:CD3E"）を color に指定する sc.pl.umap() もサポートされています。

import matplotlib.pyplot as plt

# Joint UMAP colored by cluster assignment
sc.pl.umap(mdata, color="leiden_wnn", title="WNN Leiden clusters",
           legend_loc="on data", show=False)
plt.savefig("wnn_umap_clusters.png", dpi=150, bbox_inches="tight")
plt.close()
print("Saved wnn_umap_clusters.png")

# Color by RNA gene expression on joint UMAP
sc.pl.umap(mdata, color=["rna:CD3E", "rna:CD19", "rna:CD14"],
           use_raw=False, vmax="p99", show=False, ncols=3)
plt.savefig("wnn_umap_markers.png", dpi=150, bbox_inches="tight")
plt.close()
print("Saved wnn_umap_markers.png")

# Per-modality scatter / embedding plots
mu.pl.embedding(mdata, basis="X_umap", color="leiden_wnn",
                show=False)
plt.savefig("embedding_clusters.png", dpi=150, bbox_inches="tight")
plt.close()

# Violin plot: RNA QC metrics per cluster
sc.pl.violin(mdata["rna"], keys=["n_genes_by_counts", "total_counts"],
             groupby=mdata.obs.loc[mdata["rna"].obs_names, "leiden_wnn"],
             rotation=90, show=False)
plt.savefig("rna_qc_by_cluster.png", dpi=150, bbox_inches="tight")
plt.close()
print("Saved rna_qc_by_cluster.png")

Module 6: モダリティ間解析

モダリティ間解析では、モダリティをまたいで特徴量を関連付けます。たとえば、遺伝子プロモーター付近の ATAC ピークのアクセシビリティを、その遺伝子の RNA 発現と関連付けたり、MOFA+ を適用して共通の潜在因子を見つけたりします。

# --- Peak-to-gene distance annotation ---
# Annotate ATAC peaks with nearest gene (requires genomic coordinates in var)
# mdata["atac"].var should contain columns: chrom, chromStart, chromEnd
# mu.atac.tl.rank_peaks_groups(mdata, groupby="leiden_wnn")  # differential peaks

# --- Differentially accessible peaks per cluster ---
atac = mdata["atac"]
# Transfer cluster labels from MuData obs to ATAC obs
atac.obs["cluster"] = mdata.obs.loc[atac.obs_names, "leiden_wnn"].values
sc.tl.rank_genes_groups(atac, groupby="cluster", method="wilcoxon",
                        use_raw=False)
top_peaks = sc.get.rank_genes_groups_df(atac, group="0").head(5)
print("Top differential peaks in cluster 0:")
print(top_peaks[["names", "logfoldchanges", "pvals_adj"]])

# --- MOFA+ factor analysis on MuData ---
# Requires: pip install mofapy2
try:
    mu.tl.mofa(mdata, n_factors=10, seed=42,
               use_obs="all", outfile="mofa_model.hdf5")
    # Factor scores stored in mdata.obsm["X_mofa"]
    print(f"MOFA factors: {mdata.obsm['X_mofa'].shape}")
    # Variance explained per factor per modality
    # Inspect with mdata.uns["mofa"]
except ImportError:
    print("Install mofapy2 to enable MOFA factor analysis")

一般的なワークフロー

Workflow 1: 10x Multiome（RNA + ATAC）の完全な WNN 共同クラスタリング

目標：10x Genomics Multiome のペア RNA および ATAC データを処理し、WNN の統合埋め込みを作成して細胞をクラスタリングし、各クラスターの RNA マーカー遺伝子を特定します。

import muon as mu
import scanpy as sc
import numpy as np
import anndata as ad
import pandas as pd
import matplotlib.pyplot as plt

# --- 1. Simulate 10x Multiome-like data (replace with mu.read_10x_h5()) ---
np.random.seed(0)
n_cells = 400
cell_ids = [f"AAACGAATC{i:05d}-1" for i in range(n_cells)]
# Simulate two cell types with distinct expression/accessibility
labels = np.array(["typeA"] * 200 + ["typeB"] * 200)

rna_counts = np.abs(np.random.negative_binomial(6, 0.35, (n_cells, 2000))).astype(float)
rna_counts[:200, :100] += 10   # typeA marker genes
rna_counts[200:, 100:200] += 10  # typeB marker genes

atac_counts = np.abs(np.random.binomial(1, 0.05, (n_cells, 50000))).astype(float)
atac_counts[:200, :5000] += 1   # typeA accessible peaks
atac_counts[200:, 5000:10000] += 1  # typeB accessible peaks

rna  = ad.AnnData(rna_counts,
                  obs=pd.DataFrame({"cell_type": labels}, index=cell_ids),
                  var=pd.DataFrame(index=[f"gene_{j}" for j in range(2000)]))
atac = ad.AnnData(atac_counts,
                  obs=pd.DataFrame({"cell_type": labels}, index=cell_ids),
                  var=pd.DataFrame(index=[f"peak_{j}" for j in range(50000)]))

mdata = mu.MuData({"rna": rna, "atac": atac})
print(f"Loaded: {mdata}")

# --- 2. RNA preprocessing ---
sc.pp.filter_cells(mdata["rna"], min_genes=50)
sc.pp.filter_genes(mdata["rna"], min_cells=5)
sc.pp.normalize_total(mdata["rna"], target_sum=1e4)
sc.pp.log1p(mdata["rna"])
sc.pp.highly_variable_genes(mdata["rna"], n_top_genes=500)
sc.pp.pca(mdata["rna"], n_comps=30, use_highly_variable=True)
print(f"RNA PCA: {mdata['rna'].obsm['X_pca'].shape}")

# --- 3. ATAC preprocessing: TF-IDF + LSI ---
sc.pp.filter_cells(mdata["atac"], min_genes=100)
sc.pp.filter_genes(mdata["atac"], min_cells=5)
mu.atac.pp.tfidf(mdata["atac"])
mu.atac.tl.lsi(mdata["atac"], n_comps=30)
print(f"ATAC LSI: {mdata['atac'].obsm['X_lsi'].shape}")

# --- 4. Per-modality neighbors (optional for modality-specific plots) ---
mu.pp.neighbors(mdata["rna"],  use_rep="X_pca", n_neighbors=15, key_added="neighbors")
mu.pp.neighbors(mdata["atac"], use_rep="X_lsi",  n_neighbors=15, key_added="neighbors")

# --- 5. WNN joint neighbor graph ---
mu.pp.intersect_obs(mdata)   # align obs after per-modality filtering
mu.pp.neighbors(mdata, key_added="wnn",
                use_rep={"rna": "X_pca", "atac": "X_lsi"},
                n_neighbors=15, random_state=42)

# --- 6. UMAP + Leiden clustering ---
sc.tl.umap(mdata, neighbors_key="wnn", random_state=42)
sc.tl.leiden(mdata, neighbors_key="wnn", resolution=0.5, key_added="leiden_wnn")
print(f"Clusters: {mdata.obs['leiden_wnn'].value_counts().to_dict()}")

# --- 7. Marker genes per cluster ---
sc.tl.rank_genes_groups(mdata["rna"],
                        groupby=mdata.obs.loc[mdata["rna"].obs_names, "leiden_wnn"],
                        method="wilcoxon", use_raw=False)
top_markers = sc.get.rank_genes_groups_df(mdata["rna"], group="0").head(5)
print("Top RNA markers for cluster 0:")
print(top_markers[["names", "logfoldchanges", "pvals_adj"]])

# --- 8. Visualization ---
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
sc.pl.umap(mdata, color="leiden_wnn", ax=axes[0], show=False, title="WNN clusters")
sc.pl.umap(mdata, color="rna:gene_0",  ax=axes[1], show=False, title="gene_0 expression",
           use_raw=False, vmax="p99")
plt.tight_layout()
plt.savefig("multiome_wnn_pipeline.png", dpi=150, bbox_inches="tight")
plt.close()
print("Saved multiome_wnn_pipeline.png")

Workflow 2: CITE-seq（RNA + 表面タンパク質）解析

目標：ペアの RNA および抗体由来タグ（ADT／タンパク質）カウントを含む CITE-seq データを処理します。タンパク質を CLR で正規化し、タンパク質マーカーから細胞型をアノテーションして、統合 UMAP を作成します。

import muon as mu
import scanpy as sc
import numpy as np
import anndata as ad
import pandas as pd
import matplotlib.pyplot as plt
from scipy.sparse import csr_matrix

np.random.seed(1)
n_cells = 350
cell_ids = [f"CITESEQ_{i:04d}" for i in range(n_cells)]

# Simulate RNA + protein modalities
# Three cell types: T-cells (CD3+CD4+), B-cells (CD19+CD20+), Monocytes (CD14+CD16+)
n_t, n_b, n_mono = 120, 130, 100
labels = ["Tcell"] * n_t + ["Bcell"] * n_b + ["Monocyte"] * n_mono

rna_mat = np.abs(np.random.negative_binomial(4, 0.4, (n_cells, 1500))).astype(float)
rna_mat[:n_t, :50] += 15       # T-cell RNA markers
rna_mat[n_t:n_t+n_b, 50:100] += 15   # B-cell RNA markers
rna_mat[n_t+n_b:, 100:150] += 15     # Monocyte RNA markers

# 20 surface proteins: first 5 T-cell, next 5 B-cell, next 5 Monocyte, rest baseline
prot_mat = np.abs(np.random.normal(2, 0.5, (n_cells, 20)))
prot_mat[:n_t, :5] += 8          # CD3, CD4, CD5, CD7, CD8 for T-cells
prot_mat[n_t:n_t+n_b, 5:10] += 8  # CD19, CD20, CD22, CD24, CD79 for B-cells
prot_mat[n_t+n_b:, 10:15] += 8    # CD14, CD16, CD64, CD11b, HLA-DR for Monocytes

protein_names = ["CD3", "CD4", "CD5", "CD7", "CD8",
                 "CD19", "CD20", "CD22", "CD24", "CD79a",
                 "CD14", "CD16", "CD64", "CD11b", "HLA-DR",
                 "CD25", "CD56", "CD45RA", "CD45RO", "IgG-ctrl"]

rna = ad.AnnData(csr_matrix(rna_mat),
                 obs=pd.DataFrame({"cell_type": labels}, index=cell_ids),
                 var=pd.DataFrame(index=[f"gene_{j}" for j in range(1500)]))
prot = ad.AnnData(prot_mat,
                  obs=pd.DataFrame({"cell_type": labels}, index=cell_ids),
                  var=pd.DataFrame(index=protein_names))

mdata = mu.MuData({"rna": rna, "protein": prot})
print(mdata)

# --- RNA preprocessing ---
sc.pp.normalize_total(mdata["rna"], target_sum=1e4)
sc.pp.log1p(mdata["rna"])
sc.pp.highly_variable_genes(mdata["rna"], n_top_genes=500)
sc.pp.pca(mdata["rna"], n_comps=30, use_highly_variable=True)

# --- Protein CLR normalization (Centered Log-Ratio) ---
# CLR normalizes each protein across cells: log(x / geometric_mean(x))
mu.prot.pp.clr(mdata["protein"])
print("Protein CLR normalization applied")
print(f"Protein data range: [{mdata['protein'].X.min():.2f}, {mdata['protein'].X.max():.2f}]")

# PCA on protein modality
sc.pp.pca(mdata["protein"], n_comps=min(15, prot.n_vars - 1))

# --- WNN joint embedding ---
mu.pp.neighbors(mdata, key_added="wnn",
                use_rep={"rna": "X_pca", "protein": "X_pca"},
                n_neighbors=20, random_state=0)
sc.tl.umap(mdata, neighbors_key="wnn", random_state=0)
sc.tl.leiden(mdata, neighbors_key="wnn", resolution=0.4, key_added="leiden_wnn")

# --- Protein-based cell type annotation ---
# Compute mean CLR protein per cluster
prot_df = pd.DataFrame(
    mdata["protein"].X,
    index=mdata["protein"].obs_names,
    columns=protein_names
)
prot_df["cluster"] = mdata.obs.loc[mdata["protein"].obs_names, "leiden_wnn"].values
cluster_means = prot_df.groupby("cluster").mean()
print("\nMean CLR protein per cluster:")
print(cluster_means[["CD3", "CD4", "CD19", "CD14"]].round(2))

# --- Visualization ---
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
sc.pl.umap(mdata, color="leiden_wnn",   ax=axes[0], show=False, title="WNN Leiden")
sc.pl.umap(mdata, color="protein:CD3",  ax=axes[1], show=False, title="CD3 (CLR)",
           use_raw=False, vmax="p99")
sc.pl.umap(mdata, color="protein:CD19", ax=axes[2], show=False, title="CD19 (CLR)",
           use_raw=False, vmax="p99")
plt.tight_layout()
plt.savefig("citeseq_joint_umap.png", dpi=150, bbox_inches="tight")
plt.close()
print("Saved citeseq_joint_umap.png")

# --- Dot plot: protein markers per cluster ---
sc.pl.dotplot(mdata["protein"],
              var_names=["CD3", "CD4", "CD19", "CD20", "CD14", "CD16"],
              groupby=mdata.obs.loc[mdata["protein"].obs_names, "leiden_wnn"],
              show=False)
plt.savefig("citeseq_protein_dotplot.png", dpi=150, bbox_inches="tight")
plt.close()
print("Saved citeseq_protein_dotplot.png")

主なパラメーター

パラメーター	モジュール／関数	デフォルト	範囲／オプション	効果
`n_neighbors`	`mu.pp.neighbors()`	`30`	`10`–`100`	グラフの近傍サイズ。値が小さいほど局所構造が細かくなります
`use_rep`	`mu.pp.neighbors()`	`None`（自動）	`{modality: embedding_key}` の辞書	WNN で使用する各モダリティの埋め込み
`key_added`	`mu.pp.neighbors()`	`"neighbors"`	任意の文字列	グラフを保存するキー。統合グラフには `"wnn"` を使用します
`n_comps`	`sc.pp.pca()`、`mu.atac.tl.lsi()`	`50`	`10`–`100`	次元削減後の次元数。通常は 30～50 です
`resolution`	`sc.tl.leiden()`	`1.0`	`0.1`–`2.0`	クラスタリングの粒度。値が小さいほど、クラスター数が少なく、規模が大きくなります
`scale_factor`	`mu.atac.pp.tfidf()`	`1e4`	`1e3`–`1e5`	対数変換前の TF スケーリング定数
`n_factors`	`mu.tl.mofa()`	`10`	`5`–`50`	MOFA の潜在因子数。説明分散のエルボーを参考にしてください
`target_sum`	`sc.pp.normalize_total()`	`1e4`	`1e3`–`1e6`	細胞あたりのライブラリーサイズ正規化の目標値（RNA）
`n_top_genes`	`sc.pp.highly_variable_genes()`	変動	`1000`–`5000`	PCA に使用する、高変動遺伝子の保持数

ベストプラクティス

WNN の前に obs を揃える：モダリティごとの QC フィルタリング後、モダリティ間で細胞数が異なる場合があります。WNN の前に必ず mu.pp.intersect_obs(mdata) を呼び出し、すべてのモダリティで同じ細胞セットを使用してください。
```
mu.pp.intersect_obs(mdata)
print(f"Shared cells after intersect: {mdata.n_obs}")
```
ATAC の LSI 第 1 成分を除外する：最初の LSI 成分は生物学的シグナルではなく、シーケンス深度（細胞あたりの総カウント）を捉えます。mu.atac.tl.lsi() はすべての成分を X_lsi に保存するため、WNN で X_lsi を使用する前に、第 1 成分が log_total_counts と相関するか確認してください。相関する場合は、使用する範囲を X_lsi[:, 1:] に限定します。
```
import numpy as np, pandas as pd
atac = mdata["atac"]
corr = np.corrcoef(atac.obsm["X_lsi"][:, 0],
                   np.log1p(atac.obs["total_counts"]))[0, 1]
print(f"LSI1 vs log_counts correlation: {corr:.3f}")
# If |corr| > 0.9, exclude component 1:
# atac.obsm["X_lsi"] = atac.obsm["X_lsi"][:, 1:]
```
key_added="wnn" を一貫して使用する：統合グラフには常に "wnn" という名前を付け、後続のすべての sc.tl.umap()、sc.tl.leiden()、sc.tl.paga() の呼び出しに neighbors_key="wnn" を渡してください。neighbors_key の値を混在させると、気付かないうちに誤ったグラフが使用されます。
タンパク質にはログ正規化ではなく CLR 正規化を使用する：CITE-seq の抗体由来タグ（ADT）カウントは、RNA とは異なるノイズモデルを持ちます。タンパク質ごとの CLR 正規化には mu.prot.pp.clr() を使用してください。タンパク質モダリティに sc.pp.normalize_total() と sc.pp.log1p() を適用しないでください。
正規化の前に RNA の生カウントを保存する：正規化を行う前に、RNA の生カウントを mdata["rna"].layers["counts"] に保存し、mdata["rna"].raw = mdata["rna"] を設定してください。これは、後続の差次的発現検定と scvi-tools との統合に必要です。
```
import scipy.sparse as sp
mdata["rna"].layers["counts"] = mdata["rna"].X.copy()
# Store pre-normalization snapshot
mdata["rna"].raw = mdata["rna"]
```
WNN 埋め込みの次元数を揃える：WNN に渡す RNA の PCA と ATAC の LSI 埋め込みは、同程度の次元数（例：どちらも 30 または 50 成分）にしてください。次元数が異なってもエラーにはなりませんが、次元数の少ないモダリティの重みが小さくなる場合があります。

よく使うレシピ

レシピ：WNNからクラスターごとのモダリティ重みを計算する

使用する場面：WNNグラフで、どの細胞がRNAシグナルとATACシグナルのどちらにより依存しているかを調べます。クラスター内でRNAの重みが高い場合はRNAデータがよりクリーンであることを示し、ATACの重みが高い場合はクロマチンアクセシビリティシグナルがより強いことを示します。

# WNN stores per-cell modality weights in mdata.obsm after mu.pp.neighbors()
# Key name: "{key_added}_weights" — check available keys
print([k for k in mdata.obsm.keys() if "wnn" in k.lower()])

# If weights are stored (depends on muon version):
if "wnn_weights" in mdata.obsm:
    weights_df = pd.DataFrame(
        mdata.obsm["wnn_weights"],
        index=mdata.obs_names,
        columns=["rna_weight", "atac_weight"]
    )
    weights_df["cluster"] = mdata.obs["leiden_wnn"].values
    print(weights_df.groupby("cluster").mean().round(3))
else:
    # Proxy: compare RNA vs ATAC PCA explained variance per cell
    rna_var  = np.var(mdata["rna"].obsm["X_pca"],  axis=1)
    atac_var = np.var(mdata["atac"].obsm["X_lsi"], axis=1)
    mdata.obs["rna_signal_proxy"]  = rna_var / (rna_var + atac_var)
    mdata.obs["atac_signal_proxy"] = atac_var / (rna_var + atac_var)
    print(mdata.obs[["rna_signal_proxy", "atac_signal_proxy", "leiden_wnn"]].groupby("leiden_wnn").mean().round(3))

レシピ：クラスターラベルを各モダリティのAnnDataに書き戻す

使用する場面：WNNによる統合クラスタリングの後、共有のクラスターラベルを各モダリティの.obsにコピーして、モダリティ固有の下流解析（例：クラスター内の差次的ピーク解析）に使用します。

# Copy joint cluster labels into each modality's obs
for mod_name in mdata.mod:
    mod_adata = mdata[mod_name]
    shared_cells = mod_adata.obs_names
    mod_adata.obs["leiden_wnn"] = mdata.obs.loc[shared_cells, "leiden_wnn"].values
    print(f"{mod_name}: added leiden_wnn, {mod_adata.obs['leiden_wnn'].nunique()} clusters")

# Now run modality-specific differential analysis per cluster
sc.tl.rank_genes_groups(mdata["rna"],  groupby="leiden_wnn",
                        method="wilcoxon", use_raw=False)
sc.tl.rank_genes_groups(mdata["atac"], groupby="leiden_wnn",
                        method="wilcoxon", use_raw=False)
print("Differential features computed for both modalities")

レシピ：MuDataをscvi-tools用の連結AnnDataに変換する

使用する場面：scvi-toolsのMultiVIとtotalVIでは、varでモダリティの識別情報を追跡する連結AnnDataが必要です。このレシピでは、これらのモデルに入力するデータを準備します。

# Concatenate RNA and ATAC into a single AnnData with modality column in var
rna_adata  = mdata["rna"].copy()
atac_adata = mdata["atac"].copy()

rna_adata.var["modality"]  = "Gene Expression"
atac_adata.var["modality"] = "Peaks"

# Restore raw counts for scvi-tools (requires integer counts)
# Ensure mdata["rna"].layers["counts"] and mdata["atac"].X are integer
import anndata as ad
combined = ad.concat([rna_adata, atac_adata], axis=1, merge="unique")
combined.obs = mdata.obs.loc[combined.obs_names].copy()
print(f"Combined AnnData: {combined.n_obs} cells × {combined.n_vars} features")
print(combined.var["modality"].value_counts())
# Use combined with scvi.model.MULTIVI.setup_anndata(combined, ...)

トラブルシューティング

問題	原因	解決策
MuDataの構築時に`KeyError: 'rna'`が発生する	モダリティキーが設定されていないか、想定された名前がないファイルから読み込まれている	`mdata.mod.keys()`を確認してください。10x h5ファイルを読み込むときは、モダリティに自動で`"rna"`と`"atac"`という名前が付く`mu.read_10x_h5()`を使用してください
WNNの実行中に`ValueError: obs_names mismatch`が発生する	モダリティごとに細胞が異なる方法でフィルタリングされ、obsが不一致になっている	モダリティごとのQC後に`mu.pp.intersect_obs(mdata)`を呼び出し、細胞セットを揃えてください
UMAP／Leidenが誤ったグラフを使用する	複数の近傍グラフが存在し、関数がデフォルトキー`"neighbors"`を使用している	`sc.tl.umap()`と`sc.tl.leiden()`には、常に`neighbors_key="wnn"`を明示的に渡してください
LSIの第1成分がATAC埋め込みを支配する	第1成分が生物学的特徴ではなく、シーケンス深度を捉えている	`X_lsi[:, 0]`と`log_total_counts`の相関を計算してください。
`mu.atac.pp.tfidf`で疎行列エラーが発生する	入力行列が密行列であるか、データ型が正しくない	先に次のように変換してください：`mdata["atac"].X = scipy.sparse.csr_matrix(mdata["atac"].X.astype(float))`
CLR正規化でタンパク質にNaNが生じる	すべてのタンパク質のカウントが0の細胞がある	CLRの前に`sc.pp.filter_cells(mdata["protein"], min_genes=1)`で細胞をフィルタリングしてください
`mu.tl.mofa()`でmofapy2が見つからず失敗する	mofapy2がインストールされていない	`pip install mofapy2`を実行してください。MOFAとの統合には`muon`のバージョンが0.1.5以上であることを確認してください
大規模なATACピーク行列でメモリエラーが発生する	ATAC行列（50k以上のピーク × 10k以上の細胞）がRAMの容量を超えている	`sc.pp.highly_variable_genes()`を使用して上位50kのピークを先に選択するか、チャンク単位で読み込んでください

参考資料

muonのドキュメント — 公式ドキュメント、チュートリアル、APIリファレンス
muon GitHub（scverse） — ソースコード、課題、使用例
Bredikhin et al. Genome Biology 2022 — muon／MuDataの論文
10x Genomics Multiome解析ガイド — muonを使ったMultiomeの手順別チュートリアル
WNN手法（Hao et al. Cell 2021） — Seurat v4のWeighted Nearest Neighborに関する原著論文

Markdown 対応

ログインしてコメント

基本情報

作者: jaechang-hits
リポジトリ: jaechang-hits/sciagent-skills
カテゴリ: データ・データベース
インストール数: 79
GitHub スター: 384

問題を報告