scrna-seq-pipeline

scRNA-seq パイプライン（Scanpy）

概要

このスキルでは、Scanpy と AnnData を使用して、再現可能な scRNA-seq 解析を Python で実行します。標準的なワークフローである、カウントデータの読み込み → QC フィルタリング → 正規化／対数変換 → HVG → スケーリング／PCA → 近傍／UMAP → Leiden → 差次的発現解析に加え、Harmony または scVI によるバッチ補正も扱います。デフォルトパラメーターは Scanpy の PBMC3k チュートリアルの慣例に準拠しており、ほとんどの 10X Genomics データセットで安全な開始点として使えます。

すべての中間埋め込み、クラスタラベル、DE 結果を AnnData オブジェクトに格納した、単一の .h5ad ファイルを出力します。

使用方法

依存関係をインストールします。

pip install "scanpy>=1.10" anndata leidenalg python-igraph harmonypy scvi-tools

パイプラインを最初から最後まで実行します。

import scanpy as sc
import scanpy.external as sce

sc.settings.verbosity = 2
sc.settings.set_figure_params(dpi=100, facecolor="white")

# 1. Load
adata = sc.read_10x_h5("filtered_feature_bc_matrix.h5")
# Alternatives:
# adata = sc.read_10x_mtx("filtered_feature_bc_matrix/", var_names="gene_symbols", cache=True)
# adata = sc.read_h5ad("input.h5ad")
adata.var_names_make_unique()

# 2. QC
adata.var["mt"] = adata.var_names.str.startswith("MT-")
sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], percent_top=None, log1p=False, inplace=True)
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
adata = adata[adata.obs.n_genes_by_counts < 5000, :]
adata = adata[adata.obs.pct_counts_mt < 20, :].copy()

# 3. Normalize + HVG
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
adata.raw = adata  # freeze log-normalized full gene set BEFORE HVG subset
sc.pp.highly_variable_genes(adata, n_top_genes=2000, flavor="seurat")
adata = adata[:, adata.var.highly_variable].copy()

# 4. Scale + PCA + neighbors + UMAP + Leiden
sc.pp.scale(adata, max_value=10)
sc.tl.pca(adata, n_comps=50, svd_solver="arpack")
sc.pp.neighbors(adata, n_neighbors=15, n_pcs=30)
sc.tl.umap(adata)
sc.tl.leiden(adata, resolution=0.5)

# 5. Differential expression per cluster
sc.tl.rank_genes_groups(adata, groupby="leiden", method="wilcoxon")

# 6. Save
adata.write("result.h5ad", compression="gzip")

バッチ補正 — Harmony

# Assumes adata.obs["batch"] exists and PCA is already computed.
sce.pp.harmony_integrate(adata, key="batch")  # writes adata.obsm["X_pca_harmony"]
sc.pp.neighbors(adata, n_neighbors=15, n_pcs=30, use_rep="X_pca_harmony")
sc.tl.umap(adata)
sc.tl.leiden(adata, resolution=0.5)

バッチ補正 — scVI

scVI には生カウントが必要です。正規化の前に実行するか、生カウントのレイヤーを保持してください。

import scvi

# adata_raw holds unnormalized integer counts; adata.obs["batch"] is the batch key.
scvi.model.SCVI.setup_anndata(adata_raw, batch_key="batch")
model = scvi.model.SCVI(adata_raw, n_latent=30)
model.train(max_epochs=400, early_stopping=True)
adata.obsm["X_scVI"] = model.get_latent_representation()

sc.pp.neighbors(adata, n_neighbors=15, use_rep="X_scVI")
sc.tl.umap(adata)
sc.tl.leiden(adata, resolution=0.5)

応答形式

• まず、ユーザーが必要とするコマンドまたはコードを提示し、その後に説明します
• 次の構成にします：入力の確認 → 動作するコード → 主要パラメーターの説明 → 注意点
• タスクごとに、完全に動作する例を 1 つ示します。すべての選択肢は示しません
• コードコメントは簡潔で実用的なものにします（存在理由ではなく、何をするかを記述します）
• 目安：コード 50～100 行と、簡潔な説明

主要な概念

• AnnData の構成。 adata.X は現在の発現行列（細胞 × 遺伝子）です。adata.obs には細胞ごとのメタデータ（QC 指標、クラスタラベル、バッチ）が格納されます。adata.var には遺伝子ごとのメタデータ（HVG フラグ、mt）が格納されます。埋め込みは adata.obsm（X_pca、X_umap、X_scVI）に格納されます。近傍グラフと DE 結果は adata.uns に格納されます。
• .raw のスナップショット。 adata.raw = adata は、HVG に絞り込む前の、全遺伝子を含む対数正規化行列を保存します。下流のプロット関数（sc.pl.umap(..., color="GENE")、sc.pl.rank_genes_groups）は .raw を参照するため、adata を HVG に絞り込んだ後でも、任意の遺伝子を可視化できます。
• QC 閾値はデータセットによって異なります。 ヒトの 10X PBMC のようなデータでは、n_genes ∈ [200, 5000] および pct_mt < 20 が妥当なデフォルトです。設定を確定する前に、分布を確認してください（sc.pl.violin）。ヒトでは MT-、マウスでは mt- を使用します。
• HVG の flavor は重要です。 対数正規化済みデータでは flavor="seurat"（デフォルト）を使用します。生カウントでは flavor="seurat_v3" を使用します。これは整数カウントを前提としているため、対数変換済みデータではエラーになったり、誤解を招く結果になったりします。
• 近傍パラメーター。 n_neighbors は局所構造と大域構造を制御します（10～50。値が大きいほど UMAP は滑らかになります）。n_pcs は sc.pl.pca_variance_ratio のエルボーを捉える値にします（通常 20～50）。Leiden の resolution はクラスタの粒度を制御します（0.2～1.5。値が高いほどクラスタ数が増えます）。
• バッチ補正方法の選択。 Harmony は高速で、PCA 上で動作し、軽度から中程度のバッチ効果に適しています。scVI は深層生成モデルであり、生カウントが必要です。高速化には GPU が必要ですが、より強い技術的変動やマルチモーダル設計に対応できます。

クイックリファレンス

読み込み

入力	呼び出し
10X HDF5	sc.read_10x_h5("filtered_feature_bc_matrix.h5")
10X MTX ディレクトリ	sc.read_10x_mtx("path/", var_names="gene_symbols", cache=True)
既存の AnnData	sc.read_h5ad("input.h5ad")
読み込み後は常に実行	adata.var_names_make_unique()

主要パラメーターの範囲

パラメーター	デフォルト	一般的な範囲	備考
filter_cells(min_genes=)	200	100–500	空または品質の低いドロップレットを除外
filter_genes(min_cells=)	3	3–10	発現頻度の低い遺伝子を除外
n_genes_by_counts <	5000	2500–8000	上限値はダブレットの判定に使用
pct_counts_mt <	20	5–25	ヒトでは組織によって異なる
normalize_total(target_sum=)	1e4	1e4	CP10k。標準的な設定
n_top_genes（HVG）	2000	1000–5000	HVG が増えるとシグナルとノイズの両方が増加
scale(max_value=)	10	10	極端な z スコアをクリップ
pca(n_comps=)	50	30–100	下流の n_pcs にはエルボーを基準に選択
neighbors(n_neighbors=)	15	10–50	値が大きいほど多様体が滑らかになる
neighbors(n_pcs=)	30	20–50	n_comps 以下
leiden(resolution=)	0.5	0.2–1.5	値が高いほどクラスタ数が増える

DE と確認

sc.tl.rank_genes_groups(adata, groupby="leiden", method="wilcoxon")
sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)

# Top-N markers per cluster as a DataFrame
import pandas as pd
result = adata.uns["rank_genes_groups"]
groups = result["names"].dtype.names
markers = pd.DataFrame({g: result["names"][g][:20] for g in groups})

保存／読み込み

adata.write("result.h5ad", compression="gzip")
adata = sc.read_h5ad("result.h5ad")

よくある間違い

• 誤り： adata.raw = adata を設定する前に、HVG のみに絞り込む 正しい方法： adata = adata[:, adata.var.highly_variable].copy() を実行する前に、必ず対数正規化済みの全行列に対して .raw を設定します 理由： .raw がないと、後続の解析で非 HVG 遺伝子をプロットしたり、スコアを計算したりできなくなります

• 誤り： 対数正規化データに flavor="seurat_v3" を使用する 正しい方法： データの状態に合わせて HVG の flavor を選びます。seurat_v3 は未加工の整数カウントを、seurat（デフォルト）は対数正規化データを必要とします 理由： 誤った flavor を使うと、HVG の順位が無意味になったりエラーが発生したりして、後続のすべての処理が損なわれます

• 誤り： 10X データを読み込んだ後に adata.var_names_make_unique() を実行しない 正しい方法： 読み込み直後に adata.var_names_make_unique() を呼び出します 理由： 重複した遺伝子シンボルは、インデックス指定、HVG の選択、DE 解析で、気づかないうちに問題を引き起こします

• 誤り： 正規化／対数変換済みデータで scVI を実行する 正しい方法： scVI には未加工のカウントを含む AnnData を渡し、setup_anndata で batch_key を設定します 理由： scVI はカウントを直接モデル化します。正規化済みの入力は分布に関する前提に反し、不正確な潜在空間を生成します

• 誤り： マウスのミトコンドリア遺伝子に MT- プレフィックスを使う（またはヒトに mt- を使う） 正しい方法： ヒトの遺伝子には MT-、マウスの遺伝子には mt- を使います 理由： プレフィックスが誤っていると、pct_counts_mt が常にゼロになり、QC フィルターが機能しません

• 誤り： .copy() を呼び出さずに AnnData をブール値でスライスする 正しい方法： フィルタリング後は必ず .copy() を使います：adata = adata[mask, :].copy() 理由： .copy() がない場合、結果はビューになります。その後のインプレース操作で警告が出たり、失敗したりします

参考資料

• Scanpy：Wolf et al. Genome Biol 2018、https://doi.org/10.1186/s13059-017-1382-0
• Scanpy チュートリアル：https://scanpy.readthedocs.io/en/stable/tutorials.html
• scVI：Gayoso et al. Nat Biotechnol 2022、https://doi.org/10.1038/s41587-021-01206-w
• Harmony：Korsunsky et al. Nat Methods 2019、https://doi.org/10.1038/s41592-019-0619-0