マルチオミクス解析：scRNA-seq 側の前処理 ─ 統合に向けて

📚 この記事について
前処理ブロックの RNA 側です。QC・ダブレット除去・正規化・HVG 選択を、最終的に統合ツール（MultiVI / scVI）へ生カウントで渡すことを見据えたコードで進めます。専門用語はその都度説明します。
🔙 前の記事：マルチオミクスのデータ構造 ─ AnnData から MuData / muon へ（前の記事）
🔜 次の記事：「scATAC-seq 側の前処理 ─ 統合に向けて」
📌 前提：AnnData / MuData の基礎（前の記事）。Python の基本。

シングルセル・マルチオミクス（single-cell multiomics／scRNA-seq と scATAC-seq の統合解析）は、RNA 側と ATAC 側をそれぞれ整えてから統合します。この記事では RNA 側の前処理を、統合（MultiVI / scVI）を見据えて進めます。

1. 統合に向けた RNA 前処理の全体像
2. QC ─ MAD で外れ値を除く
3. ダブレット除去（Scrublet）
4. 生カウントを退避してから正規化
5. HVG 選択 ─ 統合に向けて（seurat_v3）
6. scVI に渡す形に整える
まとめ
関連記事
参考文献

1. 統合に向けた RNA 前処理の全体像

RNA 側の前処理そのものは、単独の scRNA-seq 解析とほぼ同じです。違うのはゴールで、最終的に統合ツール（MultiVI / scVI）へ生カウントで渡すところまで持っていきます。全体の流れを図で示します。

図：統合を見据えた RNA 前処理の流れ。QC・ダブレット除去のあと、生カウントを .layers へ退避してから .X を正規化する。HVG は生カウントで選び、scVI / MultiVI には生カウントを渡す。

💡 この記事のゴール
QC → ダブレット除去 → 生カウントの退避 → 正規化（.X）→ HVG 選択（生カウント）→ scVI に渡す、という流れです。単独解析と同じ手順でも、生カウントを壊さないことだけは徹底します。基本的な QC・正規化の考え方は scRNA-seq 前処理の全体像も参考にどうぞ。

2. QC ─ MAD で外れ値を除く

QC は、低品質な細胞（死細胞・空液滴・破砕細胞など）を除く工程です。指標は主に、総カウント数・検出遺伝子数・ミトコンドリア遺伝子の割合の3つです。固定しきい値でも切れますが、データに応じてMAD（中央絶対偏差）ベースで外れ値を除くのが頑健です（Luecken & Theis, Mol. Syst. Biol., 2019）。

図：QC の3指標（総カウント・検出遺伝子数・ミト率）の分布と、MAD ベースの外れ値カットの考え方。赤で示した外れ値を除外する。ミト率は生物学的な上限でも切る。

python
import scanpy as sc
import numpy as np

# ミトコンドリア遺伝子にフラグを立てる
adata.var["mt"] = adata.var_names.str.startswith("MT-")
sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], inplace=True)

# MAD ベースの外れ値（中央値から nmads 個ぶん外れた細胞）
def is_outlier(a, metric, nmads=5):
    x = a.obs[metric]
    mad = np.median(np.abs(x - np.median(x)))
    return (x < np.median(x) - nmads * mad) | (x > np.median(x) + nmads * mad)

adata.obs["outlier"] = (
    is_outlier(adata, "log1p_total_counts")
    | is_outlier(adata, "log1p_n_genes_by_counts")
    | (adata.obs["pct_counts_mt"] > 20)
)
adata = adata[~adata.obs["outlier"]].copy()

💡 MAD ベースの外れ値とは
中央値からの散らばりを MAD（中央絶対偏差）で測り、中央値 ± nmads × MAD の外を外れ値とします（nmads は 3〜5 が目安）。分布が歪んでいても効きやすく、固定しきい値より恣意性が減ります。ミト率だけは生物学的な上限（例：20%）で切ることが多いです。

3. ダブレット除去（Scrublet）

ダブレットは、2つの細胞が1つのバーコードに混じったもので、偽の中間状態として現れ、統合やアノテーションを乱します。Scrublet で疑わしい細胞を検出して除きます（Wolock et al., Cell Syst., 2019）。

python
# ダブレット（2細胞が1バーコードに混入）を検出・除去
sc.pp.scrublet(adata)
adata = adata[~adata.obs["predicted_doublet"]].copy()

⚠️ 注意点：QC はモダリティごとに
ダブレット除去・QC は各モダリティで別々に行います。ペアデータでは、このあと両モダリティを通過した細胞にそろえます（合同 QC。データ構造の記事の mu.pp.intersect_obs を参照）。

4. 生カウントを退避してから正規化

ここが統合向け前処理の肝です。正規化で .X を書き換える前に、生カウントを .layers に退避します。MultiVI / scVI は生カウントを入力に取り、正規化・バッチ補正を内部で行う設計だからです（Lopez et al., 2018／Ashuach et al., 2023）。

python
# 正規化の前に、生カウントを退避（scVI / MultiVI 用）
adata.layers["counts"] = adata.X.copy()

# .X は可視化・マーカー確認用に正規化（生は .layers に残る）
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)

⚠️ 注意点：正規化した .X をモデルに渡さない
normalize_total / log1p は .X を上書きします。先に .layers[‘counts’] へコピーしておけば、可視化・マーカー確認は正規化した .X で、モデル入力は生カウントで、と1つのオブジェクトで両立できます。

5. HVG 選択 ─ 統合に向けて（seurat_v3）

HVG（高変動遺伝子）は、細胞間で発現が大きくばらつく遺伝子で、細胞の違いを担う情報源です。全遺伝子ではなく HVG に絞ると、ノイズが減り計算も軽くなります。統合向けには、生カウントに対して seurat_v3 の手法で選び、複数サンプルがあれば batch_key でバッチ横断に選ぶのが定石です（Stuart et al., Cell, 2019）。

図：HVG 選択。平均発現に対してばらつき（分散）が大きい遺伝子（teal）を選ぶ。これらが細胞間の違いを担うため、統合・クラスタリングの入力に使う。

python
# HVG は生カウントに対して seurat_v3 で選ぶ（バッチ横断）
sc.pp.highly_variable_genes(
    adata,
    flavor="seurat_v3",
    n_top_genes=2000,
    layer="counts",
    batch_key="batch",
)

💡 なぜ seurat_v3 か（統合向け）
seurat_v3 の手法は生カウントに対して分散を安定化して HVG を選びます。scVI / MultiVI は生カウントを入力に取るため、生カウント上で選んだ HVGと相性が良いのが理由です。log 正規化後に選ぶ古い方法とは、この点で使い分けます。

6. scVI に渡す形に整える

最後に、scVI へ渡す形へ整えます。HVG に絞り、生カウントの .layers を入力に指定し、バッチ情報を渡します。scVI が学習した潜在表現は、統合後の可視化・アノテーション・下流に使えます。

python
import scvi

# HVG に絞り、生カウントを入力として scVI を設定
adata_hvg = adata[:, adata.var["highly_variable"]].copy()
scvi.model.SCVI.setup_anndata(adata_hvg, layer="counts", batch_key="batch")

model = scvi.model.SCVI(adata_hvg)
model.train()
adata.obsm["X_scVI"] = model.get_latent_representation()

⚠️ 注意点：scale / PCA を先にやらない
scVI に渡すなら、log 正規化・スケーリング・PCA を先に済ませないでください。これらはモデルが前提とする生カウントを壊します。可視化用に正規化した .X は保持してよいですが、モデルには layer=’counts’ の生カウントを渡します。

これで RNA 側は統合の準備が整いました。次回は ATAC 側の前処理を、同じく統合を見据えて進めます。

まとめ

RNA 前処理の流れは単独解析と同じでも、ゴールは「生カウントで scVI / MultiVI に渡す」。
QC は総カウント・検出遺伝子数・ミト率で、MAD ベースの外れ値除去が頑健（ミト率は上限でも切る）。
ダブレットは Scrublet で除去。QC・ダブレットはモダリティごとに。
正規化前に .layers[‘counts’] へ生カウントを退避。HVG は生カウントに seurat_v3、batch_key でバッチ横断。
scVI には HVG × 生カウント（layer=’counts’）＋ batch を渡す。scale / PCA を先にやらない。

📖 マルチオミクスのデータ構造 ─ AnnData から MuData / muon へ（前の記事） … 土台のデータ構造
📖 scRNA-seq 前処理の全体像 … RNA 単独の前処理の詳細
📖 scvi-tools の全体像 … scVI / MultiVI の土台
📖 scATAC-seq 解析の全体像 … 次の ATAC 側
📖 AnnData のデータ構造 … データ形式の基礎

参考文献

Wolf, F. A., Angerer, P., & Theis, F. J. (2018). SCANPY: large-scale single-cell gene expression data analysis. Genome Biology, 19, 15. doi:10.1186/s13059-017-1382-0
Luecken, M. D., & Theis, F. J. (2019). Current best practices in single-cell RNA-seq analysis: a tutorial. Molecular Systems Biology, 15(6), e8746. doi:10.15252/msb.20188746
Wolock, S. L., Lopez, R., & Klein, A. M. (2019). Scrublet: Computational Identification of Cell Doublets in Single-Cell Transcriptomic Data. Cell Systems, 8(4), 281–291. doi:10.1016/j.cels.2018.11.005
Stuart, T., Butler, A., Hoffman, P., et al. (2019). Comprehensive Integration of Single-Cell Data. Cell, 177(7), 1888–1902. doi:10.1016/j.cell.2019.05.031
Lopez, R., Regier, J., Cole, M. B., Jordan, M. I., & Yosef, N. (2018). Deep generative modeling for single-cell transcriptomics. Nature Methods, 15(12), 1053–1058. doi:10.1038/s41592-018-0229-2
Ashuach, T., Gabitto, M. I., Koodli, R. V., et al. (2023). MultiVI: deep generative model for the integration of multimodal data. Nature Methods, 20(8), 1222–1231. doi:10.1038/s41592-023-01909-9