scATAC-seq(single-cell ATAC-seq/核を使う場合は snATAC-seq)の解析は、まずデータを読み込むところから始まります。ここでは入力データの種類と、ATAC ならではのデータ構造を、実際のコードとともに押さえます。
1. 入力データは2種類ある
scATAC の入力は大きく2種類です。どちらから始めるかで、後の自由度が変わります。
| 入力 | 中身 | 使いどころ |
|---|---|---|
| フラグメントファイル(fragments.tsv.gz + .tbi) | 染色体・開始・終了・バーコード・カウント | 最も生に近い。ピークを後から自由に決められる(推奨) |
| 出来合いのピーク行列(.h5) | すでに決まった 細胞×ピーク | 手早く始めたいとき。ピークは固定 |
本サイトは SnapATAC2 の推奨に従い、フラグメントファイルから始めます。理由は次のデータ構造の話と直結します。
2. ATAC の行列は「タイル」と「ピーク」を区別する
ここが RNA と最も違う点です。ATAC には特徴量(=行列の列)の作り方が2通りあります。
タイル行列は、ゲノムを固定幅(例:500bp)で機械的に区切ったマス目です。ピークをまだ決めなくても作れるので、最初の次元削減の入力に使います。列数はゲノム全体で数百万にもなります。つまりリードが検出されていないゲノム座標にも、データ上はリード数「0」を当てはめます。
ピーク行列は、フラグメントが密集した意味のある領域だけを列にしたものです。ただしピークは、クラスタが決まってからでないと精度よく検出できません(希少な細胞型のピークは、全体で見ると埋もれてしまうため)。
💡 なぜ「タイルが先、ピークが後」なのか
ピークを正確に呼ぶには、まず似た細胞をまとめる(クラスタリングする)必要があります。そこで SnapATAC2 は、①固定幅のタイルで一旦クラスタリング → ②各クラスタごとにピークを検出、という順で進めます。だから最初に扱うのはタイル行列で、いきなりピーク行列ではありません(bin法の考え方:Fang et al., Nat. Commun., 2021)。
ピークを正確に呼ぶには、まず似た細胞をまとめる(クラスタリングする)必要があります。そこで SnapATAC2 は、①固定幅のタイルで一旦クラスタリング → ②各クラスタごとにピークを検出、という順で進めます。だから最初に扱うのはタイル行列で、いきなりピーク行列ではありません(bin法の考え方:Fang et al., Nat. Commun., 2021)。
3. AnnData の中身(ATAC版)
SnapATAC2 はデータを AnnData(単一細胞データの標準形式)に格納します。RNA の細胞×遺伝子に対し、ATAC は細胞×領域です。中身の対応は次の通りです。
| 場所 | 入るもの(ATAC) |
|---|---|
| obs(行=細胞) | QC指標:n_fragment・tsse・frac_dup・frac_mito |
| var(列=特徴量) | ゲノム領域。var_names は座標(例:chr1:0-500)。RNAの遺伝子名に相当 |
| X | 細胞 × 特徴量 のスパース行列(タイル or ピーク) |
| obsm[‘fragment_paired’] | フラグメントそのものの情報 |
💡 AnnData の基礎
obs / var / X といった AnnData の基本構造は、scRNA-seq 側で詳しく解説しています。ATAC で違うのは「var が遺伝子ではなくゲノム座標」という一点が中心です。
obs / var / X といった AnnData の基本構造は、scRNA-seq 側で詳しく解説しています。ATAC で違うのは「var が遺伝子ではなくゲノム座標」という一点が中心です。
4. コード:SnapATAC2 で読み込む
フラグメントを読み込み、QC指標を計算し、タイル行列を作るまでの最小コードです。
python
import snapatac2 as snap
# 1. フラグメントファイルを読み込み → backed AnnData(ディスク常駐)
data = snap.pp.import_fragments(
"fragments.tsv.gz",
chrom_sizes=snap.genome.hg38, # 参照ゲノムに合わせる
sorted_by_barcode=False,
file="atac.h5ad", # ディスクに保存(省メモリ)
)
# 2. QC指標(TSS濃縮)を計算
snap.metrics.tsse(data, snap.genome.hg38)
# 3. タイル(固定幅ビン)行列を作成 → これが次元削減の入力
snap.pp.add_tile_matrix(data, bin_size=500)
print(data)
実行後の AnnData はこうなります。n_vars ≈ 606万が、ゲノムを500bpで区切ったタイルの数です。
text
# 出力例:
AnnData object with n_obs × n_vars = 585 × 6062095
obs: 'n_fragment', 'frac_dup', 'frac_mito', 'tsse'
uns: 'reference_sequences'
obsm: 'fragment_paired'
💡 バージョンによる関数名の違い
読み込み関数は現行版では import_fragments です。古いバージョン(〜2.7)では import_data という名前でした。エラーが出たら、まず入れている SnapATAC2 のバージョンを確認してください。
読み込み関数は現行版では import_fragments です。古いバージョン(〜2.7)では import_data という名前でした。エラーが出たら、まず入れている SnapATAC2 のバージョンを確認してください。
5. 注意点
- インデックスが要る:フラグメントファイルには対応する
.tbi(tabixインデックス)が必要です。 - 参照ゲノムを一致させる:
chrom_sizesは、データを作った参照(hg38 か hg19 か)と揃えます。ズレると座標が全部おかしくなります。 - ビンサイズの目安:細胞数が多い・ざっくりでよいなら 5000bp、細かく見たいなら 500bp。大きいほど高速で粗くなります。
- backed を使う:
file=を指定してディスク常駐にしないと、数百万列の行列がメモリに乗りきらないことがあります。
まとめ
- 入力はフラグメントファイルから始めるのが柔軟(推奨)。
- 特徴量にはタイル(固定幅・最初の埋め込み用)とピーク(クラスタ後に作る)の2種類がある。
- AnnData の var は遺伝子ではなくゲノム座標。backed でディスク常駐にして省メモリに。
関連記事
- scATAC-seq 前処理の全体像
- scATAC-seq の QC
- scATAC-seq のピークコール
- scATAC-seq の次元削減とクラスタリング ─ なぜ PCA ではないのか
- scATAC-seq 解析の全体像 ─ 何を、どの順番で行うのか
参考文献
- Fang, R., Preissl, S., Li, Y., et al. (2021). Comprehensive analysis of single cell ATAC-seq data with SnapATAC. Nature Communications, 12, 1337. doi:10.1038/s41467-021-21583-9
- Zhang, K., Zemke, N. R., Armand, E. J., & Ren, B. (2024). A fast, scalable and versatile tool for analysis of single-cell omics data. Nature Methods, 21(2), 217–227. doi:10.1038/s41592-023-02139-9


コメント