single-cell-rna-qc

シングルセル RNA-seq の品質管理

scverse のベストプラクティスに従った、シングルセル RNA-seq データ向けの自動 QC ワークフローです。

このスキルを使用するタイミング

次の場合に使用します。

• ユーザーからシングルセル RNA-seq データの品質管理または QC を依頼された場合
• 低品質な細胞をフィルタリングしたい、またはデータ品質を評価したい場合
• QC の可視化やメトリクスが必要な場合
• scverse／scanpy のベストプラクティスに従うよう求められた場合
• MAD に基づくフィルタリングや外れ値検出を依頼された場合

対応する入力形式：

• .h5ad ファイル（scanpy／Python ワークフローの AnnData 形式）
• .h5 ファイル（10X Genomics Cell Ranger の出力）

デフォルトの推奨事項：ユーザーに特別なカスタム要件がある場合や、標準以外のフィルタリングロジックを明示的に求められた場合を除き、Approach 1（完全なパイプライン）を使用してください。

Approach 1：完全な QC パイプライン（標準的なワークフローに推奨）

scverse のベストプラクティスに従った標準的な QC には、便利なスクリプト scripts/qc_analysis.py を使用してください。

python3 scripts/qc_analysis.py input.h5ad
# or for 10X Genomics .h5 files:
python3 scripts/qc_analysis.py raw_feature_bc_matrix.h5

スクリプトはファイル形式を自動的に検出し、適切な方法で読み込みます。

このアプローチを使用するタイミング：

• すべての細胞に同じ方法を適用する、しきい値を調整可能な標準的な QC ワークフロー
• 複数のデータセットのバッチ処理
• 手早く探索的な分析を行う場合
• ユーザーが「そのまま使える」ソリューションを求めている場合

要件： anndata、scanpy、scipy、matplotlib、seaborn、numpy

パラメーター：

コマンドラインパラメーターを使用して、フィルタリングのしきい値と遺伝子パターンをカスタマイズします。

• --output-dir - 出力ディレクトリ
• --mad-counts、--mad-genes、--mad-mt - counts／genes／MT% の MAD しきい値
• --mt-threshold - ミトコンドリア割合のハードカットオフ
• --min-cells - 遺伝子のフィルタリングしきい値
• --mt-pattern、--ribo-pattern、--hb-pattern - 種ごとの遺伝子名パターン

現在のデフォルト値を確認するには、--help を使用してください。

出力：

デフォルトでは、すべてのファイルが <input_basename>_qc_results/ ディレクトリに保存されます（または --output-dir で指定したディレクトリに保存されます）。

• qc_metrics_before_filtering.png - フィルタリング前の可視化
• qc_filtering_thresholds.png - MAD に基づくしきい値のオーバーレイ
• qc_metrics_after_filtering.png - フィルタリング後の品質メトリクス
• <input_basename>_filtered.h5ad - 下流解析にすぐ使える、クリーンなフィルタリング済みデータセット
• <input_basename>_with_qc.h5ad - QC アノテーションを保持した元データ

ユーザーがアクセスできるように出力をコピーする場合は、ユーザーが直接プレビューできるよう、ディレクトリ全体ではなく個々のファイルをコピーしてください。

ワークフローの手順

スクリプトは次の手順を実行します。

1. QC メトリクスを計算 - カウント深度、遺伝子検出数、ミトコンドリア／リボソーム／ヘモグロビンの含有量
2. MAD に基づくフィルタリングを適用 - counts／genes／MT% の MAD しきい値を使った、許容度の高い外れ値検出
3. 遺伝子をフィルタリング - 少数の細胞でしか検出されない遺伝子を除去
4. 可視化を生成 - しきい値を重ねた、フィルタリング前後の包括的なプロット

アプローチ 2：モジュール式の構成要素（カスタムワークフロー向け）

カスタムの解析ワークフローや標準外の要件には、scripts/qc_core.py と scripts/qc_plotting.py のモジュール式ユーティリティ関数を使用します。

# Run from scripts/ directory, or add scripts/ to sys.path if needed
import anndata as ad
from qc_core import calculate_qc_metrics, detect_outliers_mad, filter_cells
from qc_plotting import plot_qc_distributions  # Only if visualization needed

adata = ad.read_h5ad('input.h5ad')
calculate_qc_metrics(adata, inplace=True)
# ... custom analysis logic here

このアプローチを使用する場合：

• 別のワークフローが必要な場合（手順を省略する、順序を変更する、サブセットごとに異なるしきい値を適用する）
• 条件分岐が必要な場合（例：他の細胞とは異なる方法でニューロンをフィルタリングする）
• 一部のみを実行する場合（フィルタリングは行わず、指標の計算や可視化のみを行う）
• 大規模なパイプライン内の他の解析手順と統合する場合
• コマンドラインパラメーターでサポートされる範囲を超えたカスタムフィルタリング条件が必要な場合

利用可能なユーティリティ関数：

qc_core.py から（コア QC 操作）：

• calculate_qc_metrics(adata, mt_pattern, ribo_pattern, hb_pattern, inplace=True) - QC 指標を計算し、adata に注釈を付けます
• detect_outliers_mad(adata, metric, n_mads, verbose=True) - MAD ベースの外れ値検出を行い、ブール値のマスクを返します
• apply_hard_threshold(adata, metric, threshold, operator='>', verbose=True) - ハードカットオフを適用し、ブール値のマスクを返します
• filter_cells(adata, mask, inplace=False) - ブール値のマスクを適用して細胞をフィルタリングします
• filter_genes(adata, min_cells=20, min_counts=None, inplace=True) - 検出に基づいて遺伝子をフィルタリングします
• print_qc_summary(adata, label='') - 要約統計量を出力します

qc_plotting.py から（可視化）：

• plot_qc_distributions(adata, output_path, title) - 包括的な QC プロットを生成します
• plot_filtering_thresholds(adata, outlier_masks, thresholds, output_path) - フィルタリングのしきい値を可視化します
• plot_qc_after_filtering(adata, output_path) - フィルタリング後のプロットを生成します

カスタムワークフローの例：

例 1：指標の計算と可視化のみを行い、まだフィルタリングしない

adata = ad.read_h5ad('input.h5ad')
calculate_qc_metrics(adata, inplace=True)
plot_qc_distributions(adata, 'qc_before.png', title='Initial QC')
print_qc_summary(adata, label='Before filtering')

例 2：MT% フィルタリングのみを適用し、他の指標は緩やかな条件のままにする

adata = ad.read_h5ad('input.h5ad')
calculate_qc_metrics(adata, inplace=True)

# Only filter high MT% cells
high_mt = apply_hard_threshold(adata, 'pct_counts_mt', 10, operator='>')
adata_filtered = filter_cells(adata, ~high_mt)
adata_filtered.write('filtered.h5ad')

例 3：サブセットごとに異なるしきい値を適用する

adata = ad.read_h5ad('input.h5ad')
calculate_qc_metrics(adata, inplace=True)

# Apply type-specific QC (assumes cell_type metadata exists)
neurons = adata.obs['cell_type'] == 'neuron'
other_cells = ~neurons

# Neurons tolerate higher MT%, other cells use stricter threshold
neuron_qc = apply_hard_threshold(adata[neurons], 'pct_counts_mt', 15, operator='>')
other_qc = apply_hard_threshold(adata[other_cells], 'pct_counts_mt', 8, operator='>')

ベストプラクティス

1. フィルタリングは緩やかに行う - デフォルトのしきい値は、希少な集団を失わないよう、意図的に大部分の細胞を保持します
2. 可視化を確認する - フィルタリング前後のプロットを必ず確認し、フィルタリングが生物学的に妥当か確かめます
3. データセット固有の要因を考慮する - 組織によってはミトコンドリア含有量が自然に高い場合があります（例：ニューロン、心筋細胞）
4. 遺伝子注釈を確認する - ミトコンドリア遺伝子の接頭辞は生物種によって異なります（マウスでは mt-、ヒトでは MT-）
5. 必要に応じて反復する - QC パラメーターは、特定の実験や組織の種類に応じて調整が必要な場合があります

参考資料

QC の手法、パラメーターの根拠、トラブルシューティングのガイダンスについて詳しくは、references/scverse_qc_guidelines.md を参照してください。この資料には次の内容が記載されています。

• 各 QC 指標の詳細な説明と、その重要性
• MAD ベースのしきい値を用いる根拠と、固定カットオフより優れている理由
• QC の可視化（ヒストグラム、バイオリンプロット、散布図）の解釈に関するガイドライン
• 遺伝子注釈に関する生物種固有の考慮事項
• フィルタリングパラメーターを調整するタイミングと方法
• 高度な QC に関する考慮事項（環境 RNA 補正、ダブレット検出）

手法をより深く理解する必要がある場合や、QC の問題をトラブルシューティングする場合は、この資料を読み込んでください。

QC 後の次のステップ

一般的な下流解析の手順：

• 環境 RNA 補正（SoupX、CellBender）
• ダブレット検出（scDblFinder）
• 正規化（対数正規化、scran）
• 特徴量選択と次元削減
• クラスタリングと細胞種の注釈