muon-multiomics-singlecell

muon — 多模态单细胞分析

概述

muon 是一个用于多模态单细胞数据分析的 Python 框架，扩展了 AnnData 生态系统。其核心数据结构 MuData 可容纳多个 AnnData 对象（每个对象对应一种模态：RNA、ATAC、蛋白质等），并共享观测和变量轴，从而支持跨所有模态的协同操作。muon 提供特定模态的预处理流程（ATAC 的 TF-IDF 和 LSI、表面蛋白的 CLR 归一化）、用于联合降维的加权近邻（WNN）图构建，以及跨模态分析工具。它可直接与 scanpy、scvi-tools 和 MOFA+ 集成，支持完整的单细胞多组学工作流。

适用场景

• 分析 10x Genomics Multiome 数据（从同一批细胞核同时获取 RNA + ATAC）
• 处理 CITE-seq 实验（从同一批细胞获取 RNA + 表面蛋白）
• 通过 WNN 构建整合两种或更多模态信号的联合 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 文件，或按模态划分的 AnnData 对象（细胞 x 特征），且共享相同的 obs_names
• 环境：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，会让测序深度主导第一主成分；随后加权图又会因为预处理导致某个模态信息不足，而降低该模态的权重。

快速开始

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 创建与 I/O

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）会生成低维嵌入。第一个 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")

模块 5：可视化

muon 扩展了 scanpy 的绘图接口，提供了感知模态的函数。mu.pl.embedding() 可根据任意模态的特征为联合 UMAP 嵌入着色；也支持使用 color 指向带模态前缀的特征名称（例如 "rna:CD3E"）调用 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")

模块 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")

常见工作流

工作流 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")

工作流 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 要保留的高变基因数量

最佳实践

1. 在 WNN 前对齐 obs：按模态进行 QC 过滤后，各模态的细胞数可能不同。运行 WNN 前，务必调用 mu.pp.intersect_obs(mdata)，确保每个模态具有相同的细胞集合。

   mu.pp.intersect_obs(mdata)
   print(f"Shared cells after intersect: {mdata.n_obs}")
   

2. 从 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:]
   

3. 始终使用 key_added="wnn"：始终将联合图命名为 "wnn"，并在所有后续的 sc.tl.umap()、sc.tl.leiden() 和 sc.tl.paga() 调用中传入 neighbors_key="wnn"。混用不同的 neighbors_key 值会悄然导致使用错误的图。

4. 蛋白质使用 CLR 归一化，不要进行对数归一化：CITE-seq 的抗体衍生标签（ADT）计数与 RNA 具有不同的噪声模型。使用 mu.prot.pp.clr() 对每种蛋白质进行 CLR 归一化。不要对蛋白质模态使用 sc.pp.normalize_total() + sc.pp.log1p()。

5. 归一化前存储原始 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"]
   

6. 匹配 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 转换为拼接后的 AnnData，以供 scvi-tools 使用

适用场景：scvi-tools 中的 MultiVI 和 totalVI 要求使用拼接后的 AnnData，并在 var 中记录模态身份。此配方可为这些模型准备输入数据。

# 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 文件时，使用 mu.read_10x_h5()，它会自动将模态命名为 "rna" 和 "atac"
WNN 过程中出现 ValueError: obs_names mismatch	各模态经过不同方式的细胞筛选，导致 obs 不匹配	各模态完成质量控制后，调用 mu.pp.intersect_obs(mdata)，使细胞集保持一致
UMAP/Leiden 使用了错误的图	存在多个邻居图；函数使用默认键 "neighbors"	始终在 sc.tl.umap() 和 sc.tl.leiden() 中显式传入 neighbors_key="wnn"
LSI 分量 1 主导 ATAC 嵌入	第一个分量捕获的是测序深度，而非生物学信号	计算 X_lsi[:, 0] 与 log_total_counts 的相关性；若
mu.atac.pp.tfidf 引发稀疏矩阵错误	输入矩阵为稠密矩阵，或数据类型错误	先转换：mdata["atac"].X = scipy.sparse.csr_matrix(mdata["atac"].X.astype(float))
CLR 归一化导致蛋白质数据出现 NaN	某个细胞的所有蛋白质计数均为零	在进行 CLR 前，使用 sc.pp.filter_cells(mdata["protein"], min_genes=1) 筛除细胞
mu.tl.mofa() 因缺少 mofapy2 而失败	未安装 mofapy2	执行 pip install mofapy2；确保 muon 版本 ≥ 0.1.5，以支持 MOFA 集成
大型 ATAC 峰矩阵导致内存错误	ATAC 矩阵（50k+ 个峰 × 10k+ 个细胞）超出 RAM 容量	使用 sc.pp.highly_variable_genes()，先选择排名靠前的 50k 个峰；或者分块加载

相关技能

• scanpy-scrna-seq — 单模态 RNA 分析；可作为 muon 中 RNA 预处理步骤的基础
• scvi-tools-single-cell — 深度生成式多模态集成（MultiVI 用于 RNA+ATAC，totalVI 用于 CITE-seq）；需要对多个样本进行概率式批次校正时使用
• mofaplus-multi-omics — 多组学因子分析；mu.tl.mofa() 在内部调用 mofapy2，并将因子存储在 MuData 中
• anndata-data-structure — AnnData 基础知识；每个 MuData 模态都是一个标准的 AnnData 对象
• deeptools-ngs-analysis — 峰调用前，对上游 ATAC-seq BAM → bigWig 进行归一化
• macs3-peak-calling — 生成用于创建 ATAC AnnData 计数矩阵的峰 BED 文件

参考资料

• 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 中加权最近邻方法的原始论文