scanpy-scrna-seq

Scanpy によるシングルセル RNA-seq 解析

概要

Scanpy は、AnnData 形式を基盤とする、シングルセル RNA-seq データ解析向けのスケーラブルな Python ツールキットです。このスキルでは、品質管理、正規化、高変動遺伝子の選択、次元削減、クラスタリング、マーカー遺伝子の特定、細胞型のアノテーションという、標準的なエンドツーエンドのワークフローを扱います。アノテーション付きデータセットと、出版物に適した品質の可視化を作成します。

使用する場面

• シングルセル RNA-seq のカウント行列（10X Genomics、h5ad、CSV、loom）を解析する場合
• scRNA-seq データセットで品質管理フィルタリング（ミトコンドリア遺伝子の割合、遺伝子数）を行う場合
• 次元削減（PCA、UMAP、t-SNE）を実行する場合
• Leiden コミュニティ検出によって細胞クラスターを特定する場合
• クラスターごとに差次的に発現するマーカー遺伝子を特定する場合（Wilcoxon、t 検定、ロジスティック回帰）
• 既知のマーカー遺伝子パネルから細胞型をアノテーションする場合
• 軌跡推定と疑似時間解析（PAGA、拡散疑似時間）を実施する場合
• 出版物に適した品質のシングルセルプロット（ドットプロット、ヒートマップ、積み上げバイオリンプロット）を作成する場合
• 細胞型内で実験条件間の遺伝子発現を比較する場合
• エクスポートしたテーブルから汎用的な図（火山図／ヒートマップ／棒グラフ／KM）を作成するには omics-plotting SKILL を使用します。シングルセルの可視化には scanpy の sc.pl.* を使用します。
• 既存の R ワークフロー内で scRNA-seq 解析を行う場合や、Seurat 固有のアッセイタイプが必要な場合は、代わりに Seurat（R/Bioconductor）を使用します。

前提条件

• Python パッケージ：scanpy>=1.10、leidenalg、igraph、anndata
• データ要件：遺伝子発現カウント行列（細胞 x 遺伝子）。一般的な形式：10X Cell Ranger 出力（filtered_feature_bc_matrix/）、.h5ad、.h5、.csv
• 環境：Python 3.9 以降。5 万個を超える細胞を扱う場合は、16GB 以上の RAM を推奨します。

pip install "scanpy[leiden]" anndata

事前確認

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

decisions:
  - id: D1
    param: pct_counts_mt
    kind: required
    source: data
    ask: "Dying cells leak cytoplasmic RNA and read as high-mitochondrial - where should the cut sit for this tissue?"
    default: "20% (PBMC ~5%, solid tumour ~20%; propose from the observed distribution)"

  - id: D2
    param: batch_key
    kind: required
    source: data
    ask: "Which variable separates samples that were processed apart (donor, run, 10X lane), so its effect is not read as biology?"
    default: null
    skip_if: "all cells come from one sample and one processing run"

  - id: D3
    param: resolution
    kind: optional
    source: user
    ask: "How finely should cells be split - broad lineages, or subtypes within them?"
    default: "0.8 (0.1-0.3 broad lineages, 1.0-2.0 fine subtypes)"

  - id: D4
    param: n_top_genes
    kind: optional
    source: user
    ask: "How many variable genes should drive the embedding? More captures subtle states but adds noise."
    default: 2000

  - id: D5
    param: n_pcs
    kind: optional_conditional
    source: data
    depends_on: [D4]
    ask: "How many principal components carry real structure, by the variance-ratio elbow?"
    default: 40

  - id: D6
    param: min_genes, min_cells
    kind: optional_conditional
    source: data
    ask: "Where should empty droplets and undetected genes be cut off?"
    default: "cells with <200 genes, genes in <3 cells"

  - id: D6b
    param: annotationStrategy
    kind: required
    source: user
    ask: "After clustering, should clusters be named from canonical markers by hand, or handed to a reference-based annotator?"
    default: "manual marker-based within this skill; see `single-cell-annotation-guide` to choose"

  - id: D6c
    param: markerPanel
    kind: required
    source: literature
    depends_on: [D6b]
    ask: "Which canonical markers define the populations expected in this tissue and state?"
    default: null
    skip_if: "annotation delegated to a reference-based tool"

  - id: D7
    param: method (rank_genes_groups)
    kind: optional
    source: user
    ask: "Which test should rank marker genes per cluster?"
    default: "wilcoxon"

  - id: D8
    param: target_sum, flavor
    kind: optional
    source: user
    ask: "Should the standard library-size normalization and Seurat v3 highly-variable-gene selection be used?"
    default: "target_sum 1e4, flavor seurat_v3"

  - id: D9
    param: random_state
    kind: never_ask
    source: data
    reason: "Fixes the draw for reproducibility; does not change what the data supports"
    default: 0

D5 は D4 に依存します。エルボーは選択した可変遺伝子を使って計算した PCA から読み取るため、遺伝子数が変わるとその位置も変わります。D2 は単一の呼び出しではなく、パイプライン全体に影響します。バッチを指定すると、近傍グラフの構築前に補正ステージ（harmony-batch-correction または ComBat）が追加されます。

ワークフロー

ステップ 1：セットアップとデータの読み込み

scanpy の設定を行い、カウント行列を AnnData オブジェクトに読み込みます。遺伝子名が一意であることを確認してください。

import scanpy as sc
import numpy as np
import pandas as pd

# Configure settings
sc.settings.verbosity = 3  # errors=0, warnings=1, info=2, hints=3
sc.settings.set_figure_params(dpi=80, facecolor="white")

# Load data — pick the reader matching your format
adata = sc.read_10x_mtx(
    "path/to/filtered_feature_bc_matrix/",
    var_names="gene_symbols",
    cache=True,
)
# Alternatives:
# adata = sc.read_h5ad("data.h5ad")
# adata = sc.read_10x_h5("data.h5")

adata.var_names_make_unique()
print(f"Loaded: {adata.n_obs} cells x {adata.n_vars} genes")
print(f"Sparsity: {1 - adata.X.nnz / (adata.n_obs * adata.n_vars):.1%}")

ステップ 2：品質管理

ミトコンドリア遺伝子とリボソーム遺伝子にアノテーションを付け、品質管理指標を計算し、分布を可視化して、品質の低い細胞を除外します。

# Annotate gene groups
adata.var["mt"] = adata.var_names.str.startswith("MT-")    # human; use "mt-" for mouse
adata.var["ribo"] = adata.var_names.str.startswith(("RPS", "RPL"))

# Calculate QC metrics
sc.pp.calculate_qc_metrics(
    adata, qc_vars=["mt", "ribo"], percent_top=None, log1p=False, inplace=True
)

# Visualize QC distributions
sc.pl.violin(
    adata,
    ["n_genes_by_counts", "total_counts", "pct_counts_mt"],
    jitter=0.4,
    multi_panel=True,
)
sc.pl.scatter(adata, x="total_counts", y="pct_counts_mt")
sc.pl.scatter(adata, x="total_counts", y="n_genes_by_counts")

# Filter — adjust thresholds based on the QC plots above
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, :]  # remove potential doublets
adata = adata[adata.obs.pct_counts_mt < 20, :]         # remove dying cells

print(f"After QC: {adata.n_obs} cells x {adata.n_vars} genes")

ステップ 3：正規化と特徴選択

ライブラリーサイズを正規化して対数変換し、高変動遺伝子（HVG）を特定します。後で差次的発現を解析できるよう、生カウントを保存します。

# Preserve raw counts in a layer
adata.layers["counts"] = adata.X.copy()

# Normalize to 10,000 counts per cell, then log-transform
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)

# Save normalized data as .raw for visualization
adata.raw = adata

# Identify highly variable genes
sc.pp.highly_variable_genes(adata, n_top_genes=2000, flavor="seurat_v3", layer="counts")
sc.pl.highly_variable_genes(adata)

# Subset to HVGs
adata = adata[:, adata.var.highly_variable]
print(f"HVG subset: {adata.n_obs} cells x {adata.n_vars} genes")

ステップ 4：スケーリングと回帰

交絡因子の影響を回帰で取り除き、遺伝子発現値をスケーリングします。これにより、PCA の準備が整います。

# Regress out unwanted sources of variation
sc.pp.regress_out(adata, ["total_counts", "pct_counts_mt"])

# Scale to unit variance, clip extreme values
sc.pp.scale(adata, max_value=10)

print(f"Scaled matrix: mean={adata.X.mean():.4f}, std={adata.X.std():.4f}")

ステップ 5：次元削減

PCA を計算し、分散比のエルボープロットを確認してから、近傍グラフを構築し、UMAP による埋め込みを行います。

# PCA
sc.tl.pca(adata, svd_solver="arpack", n_comps=50)
sc.pl.pca_variance_ratio(adata, log=True, n_pcs=50)

# Determine n_pcs from elbow plot (typically 30-50)
n_pcs = 40

# Compute k-nearest neighbor graph
sc.pp.neighbors(adata, n_neighbors=15, n_pcs=n_pcs)

# UMAP embedding
sc.tl.umap(adata)
sc.pl.umap(adata, color=["n_genes_by_counts", "total_counts", "pct_counts_mt"])

print(f"UMAP shape: {adata.obsm['X_umap'].shape}")

ステップ 6：クラスタリング

複数の解像度で Leiden クラスタリングを実行し、可視化結果を確認して最適な粒度を選択します。

# Test multiple resolutions
for res in [0.3, 0.5, 0.8, 1.0, 1.2]:
    sc.tl.leiden(adata, resolution=res, key_added=f"leiden_res{res}", flavor="igraph", n_iterations=2)

# Compare resolutions side-by-side
sc.pl.umap(
    adata,
    color=["leiden_res0.3", "leiden_res0.5", "leiden_res0.8", "leiden_res1.0"],
    ncols=2,
    legend_loc="on data",
)

# Pick final resolution based on biological plausibility
adata.obs["leiden"] = adata.obs["leiden_res0.5"]
n_clusters = adata.obs["leiden"].nunique()
print(f"Selected resolution=0.5: {n_clusters} clusters")

ステップ 7：マーカー遺伝子の特定

生カウントに対する Wilcoxon 順位和検定を使用して、クラスターごとに差次的に発現する遺伝子を特定します。

# Rank genes per cluster
sc.tl.rank_genes_groups(adata, groupby="leiden", method="wilcoxon", use_raw=True)

# Visualization
sc.pl.rank_genes_groups(adata, n_genes=20, sharey=False)
sc.pl.rank_genes_groups_dotplot(adata, n_genes=5)
sc.pl.rank_genes_groups_heatmap(adata, n_genes=10, use_raw=True, show_gene_labels=True)

# Extract as DataFrame
markers_df = sc.get.rank_genes_groups_df(adata, group=None)
top_markers = markers_df[markers_df["pvals_adj"] < 0.05].groupby("group").head(10)
print(f"Significant markers: {len(top_markers)} genes across {top_markers['group'].nunique()} clusters")
print(top_markers[["group", "names", "logfoldchanges", "pvals_adj"]].head(15))

ステップ 8：細胞型のアノテーションとエクスポート

代表的なマーカー遺伝子に基づいて細胞型を割り当て、可視化し、すべての結果を保存します。

# Define canonical markers (example: human PBMC)
marker_genes = {
    "CD4 T cells":      ["CD3D", "CD3E", "IL7R"],
    "CD8 T cells":      ["CD8A", "CD8B", "GZMK"],
    "B cells":          ["MS4A1", "CD79A", "CD79B"],
    "NK cells":         ["NKG7", "GNLY", "KLRD1"],
    "CD14+ Monocytes":  ["CD14", "LYZ", "S100A9"],
    "FCGR3A+ Monocytes":["FCGR3A", "MS4A7"],
    "Dendritic cells":  ["FCER1A", "CST3"],
    "Platelets":        ["PPBP", "PF4"],
}

# Dot plot — rows = clusters, columns = markers
sc.pl.dotplot(adata, var_names=marker_genes, groupby="leiden", use_raw=True)

# Map clusters → cell types (adjust based on your dot plot)
cluster_to_celltype = {
    "0": "CD4 T cells",
    "1": "CD14+ Monocytes",
    "2": "B cells",
    "3": "CD8 T cells",
    "4": "NK cells",
    "5": "FCGR3A+ Monocytes",
    "6": "Dendritic cells",
    "7": "Platelets",
}
adata.obs["cell_type"] = adata.obs["leiden"].map(cluster_to_celltype).fillna("Unknown")
sc.pl.umap(adata, color="cell_type", legend_loc="on data", frameon=False)

# Save results
adata.write("results/annotated_data.h5ad", compression="gzip")
adata.obs.to_csv("results/cell_metadata.csv")
markers_df.to_csv("results/marker_genes.csv", index=False)
print(f"Saved: {adata.n_obs} cells with {adata.obs['cell_type'].nunique()} cell types")

主要なパラメーター

パラメーター	デフォルト	範囲／選択肢	効果
min_genes (filter_cells)	200	100-500	細胞ごとに検出される遺伝子の最小数です。値を小さくすると、より多くの細胞が残ります
min_cells (filter_genes)	3	3-10	遺伝子を発現する細胞の最小数です。値を大きくすると、より厳しい条件になります
pct_counts_mt threshold	20	5-30	ミトコンドリア遺伝子の最大割合です。組織によって異なります（PBMCでは5%、腫瘍では20%）
n_top_genes (HVG)	2000	1000-5000	高度可変遺伝子の数です。数を増やすと、微細な変動をより多く捉えられます
flavor (HVG)	"seurat_v3"	"seurat", "seurat_v3", "cell_ranger"	HVG検出アルゴリズムです
n_pcs (neighbors)	40	20-50	近傍グラフに使用するPCの数です。エルボープロットを確認してください
n_neighbors	15	5-50	k-NNグラフのkです。値を小さくすると、局所構造が強調されます
resolution (leiden)	0.5	0.1-2.0	クラスタリングの粒度です。値を大きくすると、クラスタ数が増えます
method (rank_genes)	"wilcoxon"	"wilcoxon", "t-test", "logreg"	DE検定です。論文ではWilcoxonが推奨されます
target_sum (normalize)	1e4	1e4-1e5	細胞あたりのライブラリーサイズの目標値です

よく使われるレシピ

レシピ：ComBatによるバッチ補正

使用する場面：複数のサンプル／バッチがあり、UMAP上でバッチ効果が明確に見られる場合です。

# Requires 'batch' column in adata.obs
sc.pp.combat(adata, key="batch")
# Re-run PCA, neighbors, UMAP, clustering after correction
sc.tl.pca(adata, svd_solver="arpack")
sc.pp.neighbors(adata, n_pcs=40)
sc.tl.umap(adata)
sc.pl.umap(adata, color=["batch", "leiden"])

レシピ：PAGA + Diffusion Pseudotimeによる軌跡推定

使用する場面：細胞が個別のクラスタではなく、発生上の連続体を形成している場合です。

# PAGA graph abstraction
sc.tl.paga(adata, groups="leiden")
sc.pl.paga(adata, color="leiden", threshold=0.03)

# Reinitialize UMAP with PAGA layout
sc.tl.umap(adata, init_pos="paga")

# Diffusion pseudotime — set root cell
adata.uns["iroot"] = np.flatnonzero(adata.obs["leiden"] == "0")[0]
sc.tl.diffmap(adata)
sc.tl.dpt(adata)
sc.pl.umap(adata, color=["dpt_pseudotime", "leiden"])

レシピ：論文掲載品質の図

使用する場面：論文やプレゼンテーション用の図を作成する場合です。

sc.settings.set_figure_params(dpi=300, frameon=False, figsize=(5, 5))

# UMAP with custom palette
sc.pl.umap(
    adata, color="cell_type",
    palette="Set2",
    legend_loc="on data",
    legend_fontsize=10,
    legend_fontoutline=2,
    title="",
    save="_celltype_publication.pdf",
)

# Stacked violin plot of top markers
flat_markers = [g for genes in marker_genes.values() for g in genes[:2]]
sc.pl.stacked_violin(
    adata, var_names=flat_markers, groupby="cell_type",
    use_raw=True, swap_axes=True,
    save="_markers_violin.pdf",
)

レシピ：条件間の差次的発現

使用する場面：特定の細胞種内で、処理群と対照群を比較する場合です。

# Subset to cell type of interest
adata_t = adata[adata.obs["cell_type"] == "CD4 T cells"].copy()

# DE between conditions
sc.tl.rank_genes_groups(adata_t, groupby="condition", groups=["treated"], reference="control", method="wilcoxon", use_raw=True)
sc.pl.rank_genes_groups_volcano(adata_t, groups=["treated"])

de_results = sc.get.rank_genes_groups_df(adata_t, group="treated")
sig_genes = de_results[de_results["pvals_adj"] < 0.05]
print(f"Significant DE genes: {len(sig_genes)}")

予想される出力

• results/annotated_data.h5ad — 埋め込み（PCA、UMAP）、クラスタラベル、細胞種の注釈、raw countsレイヤーを含む完全なAnnData
• results/cell_metadata.csv — 細胞ごとのテーブル：バーコード、クラスタ、細胞種、QC指標（n_genes、total_counts、pct_mt）
• results/marker_genes.csv — クラスタごとのDE遺伝子：遺伝子名、対数フォールド変化、調整済みp値
• QC図：バイオリンプロット（n_genes、total_counts、pct_mt）、散布図
• 埋め込み図：クラスタ、細胞種、QC指標、擬似時間で色分けしたUMAP
• マーカー図：ドットプロット、ヒートマップ、標準マーカーの積み上げバイオリンプロット

トラブルシューティング

問題	原因	解決策
ModuleNotFoundError: leidenalg	Leidenパッケージがインストールされていません	pip install leidenalg igraph
PCA中にMemoryErrorが発生する	密行列がRAMの容量を超えています	sc.pp.pca(adata, chunked=True, chunk_size=20000)を使用します
UMAPが単一の塊を示す	フィルタリングが厳しすぎるか、HVGが少なすぎます	QCのしきい値を緩和し、n_top_genesを3000-5000に増やします
クラスタが多すぎる／少なすぎる	Resolutionが適切ではありません	resolutionを0.1から2.0まで変化させ、UMAPを比較します
KeyError: 'MT-'が見つからない	種に対応するプレフィックスが誤っています	マウスではmt-（小文字）、ヒトではMT-を使用します
バッチ効果がUMAPを支配している	技術的変動が補正されていません	ComBatのレシピを適用するか、統合にHarmony／scVIを使用します
adata.raw is Noneエラー	サブセット化の前にrawを保存し忘れています	正規化の後、HVGでサブセット化する前にadata.raw = adataを設定します
マーカー遺伝子が特異的でない	Resolutionが低すぎて細胞種が統合されています	Resolutionを上げるか、logistic_regressionメソッドを使用します
regress_outが遅い	細胞数が多すぎます	regress_outを省略し、sc.pp.scale()のみを使用します（多くの場合、これで十分です）
rank_genes_groupsでValueErrorが発生する	細胞数が少なすぎるクラスタがあります	DEの前に10細胞未満のクラスタを削除します：adata = adata[adata.obs.groupby('leiden').filter(lambda x: len(x)>=10).index]

同梱リソース

このスキルには、詳細を確認するための参照ファイルが含まれています。メインのSKILL.mdでさらに詳細が必要な場合に、必要に応じて参照してください。

references/api_reference.md

モジュール別（sc.pp.、sc.tl.、sc.pl.*）に整理した、すべてのscanpy関数のクイック参照表です。完全なシグネチャ、一般的なパラメーター、AnnData構造のリファレンスを掲載しています。次の場合に使用します：特定の操作に使う正確な関数名やパラメーターが必要な場合。

references/plotting_guide.md

包括的な可視化ガイドです。QCプロット、UMAPのスタイリング、マーカー遺伝子のプロット（ドットプロット、ヒートマップ、積み上げバイオリンプロット）、軌跡プロット、複数パネルの図、論文向けのカスタマイズ、カラーパレットの推奨事項を扱います。次の場合に使用します：論文やプレゼンテーション用の図を作成する場合。

references/standard_workflow.md

パイプラインの各段階について、判断ポイント、パラメーターの根拠、代替手法（MADに基づくフィルタリング、scran正規化、自動アノテーション）を含む、詳細なステップごとのリファレンスです。次の場合に使用します：解析を一から行い、上記のワークフローセクションより詳しいガイダンスが必要な場合。

参考資料

• Scanpy documentation — 公式APIリファレンスとチュートリアル（BSD-3）
• Scanpy PBMC3k tutorial — 標準的な解説（BSD-3）
• Galaxy Training: Clustering 3K PBMCs with Scanpy — ステップごとのチュートリアル（CC-BY）
• Luecken & Theis (2019) — 「Current best practices in single-cell RNA-seq analysis: a tutorial」、Mol Syst Biol
• Heumos et al. (2023) — 「Best practices for single-cell analysis across modalities」、Nat Rev Genet
• scverse ecosystem — 関連ツール：anndata、scvi-tools、squidpy、cellrank