scRNA-seq解析：データの読み込みと品質管理（QC）

🎯 この記事の結論（先に3行）

• やること：装置が出したデータを読み込み、「使えない細胞」を統計的に見つけて除く
• 迷ったら：フィルタは固定閾値より MAD（データに決めさせる方法） が再現性◎
• 出口：きれいな AnnData → 次は統合（Harmony / scVI）

この記事のゴール

装置から出てきたデータを読み込み、「使えない細胞」を除いて、きれいなデータ（AnnData：scanpy が使うデータ形式）を作るところまでです。3つの指標で1細胞ずつチェックし、明らかに使えないもの（壊れている・空っぽ・2細胞が混ざっている）を除きます。この工程を省くと、後の解析で「存在しない細胞集団」が生まれたり、本物の違いが隠れたりします。

進め方はこの5ステップです。

1. データを読み込む
2. QC指標（3つ）を計算する
3. 指標をグラフにして、データのクセを見る
4. 使えない細胞を除く（2つのやり方を比較）
5. ノイズになる遺伝子を除き、サンプルを1つにまとめる

0. そもそもデータはどんな形？

Cell Ranger という前処理ソフトを通すと、「細胞 × 遺伝子」の巨大な表（カウント行列）が得られます。各マスは「その細胞で、その遺伝子が何分子検出されたか」という数字です。

ここで知っておきたい言葉が2つあります。

• バーコード（barcode）：液滴ごとに付く固有の識別子（ID）。同じバーコードを持つリードが「1細胞分」としてまとめられます。ただし前述のとおり、バーコード1つが生細胞1個に対応するとは限りません（空の液滴にもバーコードは付く＝空液滴、1つの液滴に2細胞が入ることもある＝ダブレット）。
• UMI（unique molecular identifier）：mRNA を1分子ずつ重複なく数えるための「分子の通し番号」です。これにより「その遺伝子が何分子あったか」を正確に数えられます。

つまり QC とは、この「バーコードと実際の細胞のズレ」を統計的に取り除く作業だと考えると分かりやすいです。

1. 準備：ライブラリの読み込み

まずは使う道具（ライブラリ）を読み込みます。下のコードはそのまま実行すればOKです。

pythonコピーimport numpy as np
import scanpy as sc
import anndata as ad
from scipy.stats import median_abs_deviation
import matplotlib.pyplot as plt

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

2. データを読み込む

2-1. まずは1サンプルを読む

Cell Ranger の出力フォルダ（matrix.mtx.gz・barcodes.tsv.gz・features.tsv.gz の3点が入ったフォルダ）を指定するだけです。

pythonコピーadata = sc.read_10x_mtx(
    "path/to/filtered_feature_bc_matrix/",   # ← 自分の出力フォルダに置き換える
    var_names="gene_symbols",                # 遺伝子シンボルを名前に
    cache=True,                              # 2回目以降の読み込みを高速化
)
print(adata)

2-2. 複数サンプルをまとめて扱う（テンプレート）

実験ではふつう複数サンプル（例：対照群×2、処理群×2）を扱います。サンプルの定義を1か所にまとめておくと、後で増減しても1行直すだけで済みます。

pythonコピー# (サンプルのラベル, 10x 出力ディレクトリ, 群ラベル)
SAMPLES = [
    ("Ctrl1", "path/to/ctrl1/filtered_feature_bc_matrix/", "Ctrl"),
    ("Ctrl2", "path/to/ctrl2/filtered_feature_bc_matrix/", "Ctrl"),
    ("Treat1", "path/to/treat1/filtered_feature_bc_matrix/", "Treat"),
    ("Treat2", "path/to/treat2/filtered_feature_bc_matrix/", "Treat"),
]

adatas = {}
for label, path, group in SAMPLES:
    a = sc.read_10x_mtx(path, var_names="gene_symbols", cache=True)
    a.var_names_make_unique()        # 重複する遺伝子名をユニーク化
    a.obs["group"] = group           # 群ラベルを付与
    adatas[label] = a
    print(f"[{label}] {a.n_obs} cells × {a.n_vars} genes")

pythonコピーimport scvelo as scv
ldata = scv.read("path/to/sample.loom", cache=True)
scv.utils.clean_obs_names(adata)   # バーコード名を正規化（形式の不一致対策）
scv.utils.clean_obs_names(ldata)
adata = scv.utils.merge(adata, ldata)   # spliced / unspliced が layers に入る

3. QCで見る3つの指標（何を見るか）

QC では、細胞ごとに次の3つの数字を計算します。それぞれが何を意味し、外れているとき何を疑うかを押さえておくと、数字を機械的に当てはめずに済みます。

4. 指標を計算して、まず「目で見る」

下のコードで3つの指標を一度に計算します。いきなり閾値を決めず、まずグラフで分布を眺めるのがコツです。

pythonコピー# 例として1サンプルで進める（実運用ではサンプルごとにループ）
adata = adatas["Ctrl1"]

# ミトコンドリア遺伝子に印をつける（マウスは "mt-"、ヒトは "MT-"）
adata.var["mt"] = adata.var_names.str.startswith("mt-")

# QC 指標を計算
sc.pp.calculate_qc_metrics(
    adata,
    qc_vars=["mt"],
    percent_top=[20],   # 発現上位20遺伝子が占める割合 → 「中身の偏り」の指標
    log1p=True,
    inplace=True,
)

計算できたら、分布をグラフにします。

pythonコピー# バイオリンプロット：3指標の分布を確認
sc.pl.violin(
    adata,
    ["n_genes_by_counts", "total_counts", "pct_counts_mt"],
    jitter=0.4, multi_panel=True,
)

# 散布図：「総カウント vs 遺伝子数」を、色でミトコンドリア%を表示
sc.pl.scatter(adata, x="total_counts", y="n_genes_by_counts", color="pct_counts_mt")

5. 使えない細胞を除く：2つのやり方（固定閾値 vs MAD）

ここがこの記事の核心です。フィルタの決め方には2つの方針があり、違いを知っておくと再現性が上がります。

5-1. 固定閾値（手軽だけど根拠が弱い）

「mt% > 5% は除く」「遺伝子数 200〜6000 だけ残す」のように数値を手で決める方法です。

pythonコピー# 固定閾値の例（あくまで例。データごとに分布を見て決める必要がある）
adata_hard = adata[
    (adata.obs.n_genes_by_counts > 200)
    & (adata.obs.n_genes_by_counts < 6000)
    & (adata.obs.pct_counts_mt < 5)
].copy()
print(f"固定閾値: {adata.n_obs} -> {adata_hard.n_obs} cells")

5-2. MAD ベース（データに決めさせる・推奨）★

そこで今のおすすめは、分布から外れ値を自動で見つける MAD（中央値絶対偏差、median absolute deviation）法です（Heumos et al., 2023）。考え方は「真ん中（中央値）から大きく外れた細胞を外れ値とみなす」だけです。

pythonコピーdef is_outlier(adata, metric: str, nmads: int = 5):
    """指定した QC 指標について、中央値 ± nmads×MAD の外を外れ値とする。"""
    M = adata.obs[metric]
    med = np.median(M)
    mad = median_abs_deviation(M)
    return (M < med - nmads * mad) | (med + nmads * mad < M)

# 一般的な指標は緩め（5 MAD）、mt% は影響が大きいので厳しめ（3 MAD）に
adata.obs["outlier"] = (
    is_outlier(adata, "log1p_total_counts", 5)
    | is_outlier(adata, "log1p_n_genes_by_counts", 5)
    | is_outlier(adata, "pct_counts_in_top_20_genes", 5)
)
adata.obs["mt_outlier"] = is_outlier(adata, "pct_counts_mt", 3) | (adata.obs.pct_counts_mt > 8)

print("品質外れ値:", int(adata.obs.outlier.sum()))
print("mt 外れ値:", int(adata.obs.mt_outlier.sum()))

# 外れ値を除く
n_before = adata.n_obs
adata = adata[(~adata.obs.outlier) & (~adata.obs.mt_outlier)].copy()
print(f"MAD ベース: {n_before} -> {adata.n_obs} cells")

ここで書いた is_outlier は、NumPy と scipy だけで組まれた数行の関数です。中央値を取り、MAD を計算し、その外側を True で返しているだけで、特別なことはしていません。ただ、この数行が読めるかどうかで、QC の自由度がはっきり変わります。読めれば「どの指標に外れ値判定をかけるか」「上側だけ切るか両側を切るか」を自分で設計できますが、読めないと nmads の数字を上下させて結果を眺めるところで手が止まります。

本記事で「QC のパラメータを何度か変えて計算し直すのが近道」と書いたのは、そのとき何を変えているのかを判断できることが前提です。『改訂 独習Pythonバイオ情報解析』は、この adata.obs への列追加や条件による行の絞り込みにあたる操作を、pandas の基本として一通り扱っています。Cell Ranger の出力ファイルを開いて中身を確かめる場面で効いてくる、ファイルの読み書きや正規表現の章もあります。こういった本で体系的に解析を学ぶことも大切です。

改訂 独習Pythonバイオ情報解析（実験医学別冊）

羊土社・2025年1月／7,150円（税込）。Miniconda での環境構築、Jupyter Notebook、Python の基本文法、NumPy・pandas による表データの操作を順に押さえたうえで、Scanpy を用いたシングルセル解析へ進む構成。文字列処理やファイルの読み書きの章もある。

※上記リンクはアフィリエイト広告（PR）です

6. 遺伝子のフィルタ＋サンプルをまとめる

細胞のフィルタが終わったら、ごく少数の細胞でしか出てこない遺伝子を除きます（ほぼノイズで、計算も重くなるため）。

pythonコピーsc.pp.filter_genes(adata, min_cells=3)   # 3細胞未満でしか出ない遺伝子を除外
print("遺伝子フィルタ後:", adata.shape)

最後に、QC 済みの各サンプルを1つの AnnData にまとめます。

pythonコピー# 各サンプルに QC を適用済みとして、まとめて結合する
combined = ad.concat(
    adatas,                # {label: AnnData} の辞書
    label="sample",        # どのサンプル由来かを obs["sample"] に記録
    join="outer",          # 遺伝子セットの和集合を取る
    index_unique="-",      # バーコード衝突を防ぐためサンプル名を付加
)
print(combined)

7. QC は「除く」だけじゃない（発展）

ここまでは「明らかに使えない細胞を除く」基本の QC でした。より厳密にやるなら、次の2つも検討します（本シリーズでは別記事で解説）。

⚠️ つまずきポイントまとめ

• mt の接頭辞を間違える（ヒト MT-／マウス mt-）→ mt% が全ゼロに。最初に確認。
• いきなり閾値を決める→ まず分布をグラフで見てから。
• 複数サンプルを混ぜて一括フィルタ→ サンプル差を巻き込む。MAD は各サンプルの分布に基づくのが基本。
• 除去しすぎ→ 本物の希少細胞型まで消えることがある。除去後にクラスタが不自然に減っていないか点検。
• QCで完璧を目指しすぎる→ 後段のアノテーションが通ればOK。通らなければパラメータを振って再試行。

まとめ

• scRNA-seq の「細胞」はノイズ（空液滴・ダブレット・死細胞）を含むので、QC が前処理の入口。
• 主役は3指標：total_counts・n_genes_by_counts・pct_counts_mt。意味を理解してから閾値を決める。
• フィルタは 固定閾値（数値で一律）vs MAD（分布から自動判定）。再現性重視なら MAD が推奨。
• 仕上げに遺伝子フィルタと anndata.concat でまとめる。

次のステップ：サンプル統合とバッチ補正 — CPU で完結する古典手法 Harmony と、GPU 必須の AI 手法 scVI を、それぞれ別記事で解説します。どちらか一方を選びます。

関連記事

• 📖 scRNA-seqの前処理 入門 ─ 全体像と進め方 … このシリーズの地図（ハブ）
• 📖 scRNA-seq解析の土台：AnnDataのデータ構造 … 本記事で作る AnnData の基礎
• 📖 サンプル統合：Harmony（古典・CPU） ／ 📖 scVI（AI・GPU） … 次のステップ
• 📖 ダブレット検出：scrublet（古典） ／ 📖 SOLO（AI） … 発展的 QC

参考文献

• Heumos, L., Schaar, A. C., Lance, C., et al. (2023). Best practices for single-cell analysis across modalities. Nature Reviews Genetics, 24, 550–572. doi:10.1038/s41576-023-00586-w
• 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
• 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
• Single-cell best practices（オンライン書籍）: https://www.sc-best-practices.org