【シングルセル】Signacを用いたsingle-cell ATAC-Seq 解析のやり方【Seurat】

Rを用いたsingle-cell ATAC-Seq data解析の方法をご存じでしょうか。single-cell ATAC-Seq data解析を行うと、各細胞ごとにクロマチンアクセシビリティー状態の分布を解析することが可能です。とても便利な手法なのでぜひチャレンジしてみましょう。

Tested Environment

windows 10, R 4.4.0, Seurat 4.3.0, EnsDb.Hsapiens.v75 2.99.0, tidyverse 2.0.0, hdf5r 1.3.11, biovizBase 1.52.0, Signac 1.13

公共データを用いたRNA-seq解析に関する初心者向け技術書を販売中

画面キャプチャをふんだんに掲載したわかりやすい解説！
自宅PCからドライ研究が始められます

技術書ページへ

画像に alt 属性が指定されていません。ファイル名: -2023-05-13-19.57.28.png

公共データを用いたSingle Cell RNA-seq解析に関する初心者向け技術書を販売中

プログラミング初心者でも始められるわかりやすい解説！
RとSeuratで始めるSingle Cell RNA-seq解析！

技術書ページへ

公共データを用いたシングルセルダイナミクス解析 ~Trajectory解析、RNA Velocity解析、空間トランスクリプトーム解析のやり方~

公共データを用いたシングルセルダイナミクス解析に関する初心者向け技術書を販売中

シングルセルデータの高度な解析であるTrajectory解析、RNA Velocity解析、空間トランスクリプトーム解析が環境構築方法から詳しく解説されています！

技術書ページへ

single-cell ATAC-Seq data解析とは？

single-cell ATAC-Seq data解析とは、single-cell Assay for Transposase-accessible Chromatin with Sequencingの略で、オープンクロマチン領域の配列を細胞ごとに解析する手法のことです。以下でもう少し詳しく説明します.

ATAC-Seqとは

ATAC-Seqでは、オープンクロマチン領域のシーケンシングを行うことで、遺伝子発現の調節に関係がある可能性がある領域の塩基配列を得ることができます。ATAC-Seqの原理について図も含め説明します。クロマチンがオープンであるとき、DNAは核小体や他の構造体から解かれ、転写因子やその他の調節タンパク質が結合できる状態になります。これにより、遺伝子発現の調節が可能になります。ATAC-Seqでは、Tn5トランスポザーゼにシークエンス用タグ配列を挿入した複合体を用いることで、オープンクロマチン領域をシーケンシングします。

Diagenode社のWebサイト（https://diagenode.co.jp/techs-used-to-study-epigenetics/atac-seq）より引用

single-cell ATAC-Seqとは

single-cell ATAC-Seq data解析とは、各細胞ごとに上述したATAC-Seqを行うことにより、細胞ごとのクロマチンアクセシビリティー状態の分布を解析することが可能です。

sc ATAC-Seq data解析の実装方法

以上述べたことを基に、こちらのGithub（https://github.com/kpatel427/YouTubeTutorials/blob/main/scATACSeq_standard_workflow.R）を参考に、single-cell ATAC-Seq data解析を行います。以下に実装例の全体を示します。

library(Signac)
library(Seurat)
library(EnsDb.Hsapiens.v75)
library(tidyverse)
library(hdf5r)
library(biovizBase)

frag.file <- read.delim('atac_v1_pbmc_10k_fragments.tsv.gz', header = F, nrows = 10)
counts <- Read10X_h5('atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5')

chrom_assay <- CreateChromatinAssay(
  counts = counts,
  sep = c(":", "-"),
  fragments = "atac_v1_pbmc_10k_fragments.tsv.gz",
  min.cells = 10,
  min.features = 200
)

metadata <- read.csv(file = 'atac_v1_pbmc_10k_singlecell.csv', header = T, row.names = 1)

pbmc <- CreateSeuratObject(
  counts = chrom_assay,
  meta.data = metadata,
  assay = 'ATAC'
)

pbmc@assays$ATAC@annotation
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v75)
seqlevels(annotations) <- paste0('chr', seqlevels(annotations))

Annotation(pbmc) <- annotations
pbmc@assays$ATAC@annotation

pbmc <- NucleosomeSignal(pbmc)

pbmc <- TSSEnrichment(object = pbmc, fast = FALSE)

pbmc$blacklist_ratio <- pbmc$blacklist_region_fragments / pbmc$peak_region_fragments
pbmc$pct_reads_in_peaks <- pbmc$peak_region_fragments / pbmc$passed_filters * 100

a1 <- DensityScatter(pbmc, x = 'nCount_ATAC', y = 'TSS.enrichment', log_x = TRUE, quantiles = TRUE)
a2 <- DensityScatter(pbmc, x = 'nucleosome_signal', y = 'TSS.enrichment', log_x = TRUE, quantiles = TRUE)

a1 | a2

VlnPlot(object = pbmc, 
        features = c('nCount_ATAC', 'nFeature_ATAC', 'TSS.enrichment', 'nucleosome_signal', 'blacklist_ratio', 'pct_reads_in_peaks'),
        pt.size = 0.1,
        ncol = 6)

pbmc <- subset(x = pbmc,
               subset = nCount_ATAC > 3000 &
                 nCount_ATAC < 30000 &
                 pct_reads_in_peaks > 15 & 
                 blacklist_ratio < 0.05 &
                 nucleosome_signal < 4 &
                 TSS.enrichment > 3)

pbmc <- RunTFIDF(pbmc) # normalization
pbmc <- FindTopFeatures(pbmc, min.cutoff = 'q0') # selecting top features
pbmc <- RunSVD(pbmc) # dimensionality reduction
DepthCor(pbmc)

pbmc <- RunUMAP(object = pbmc, reduction = 'lsi', dims = 2:30)
pbmc <- FindNeighbors(object = pbmc, reduction = 'lsi', dims = 2:30)
pbmc <- FindClusters(object = pbmc, algorithm = 3)
DimPlot(object = pbmc, label = TRUE) + NoLegend()

実行結果

最終的に以下のようなUMAPによる図が出力として得られます。

コードの解説

上に書いたソースコードの解説をしていきます。

library(Signac)
library(Seurat)
library(EnsDb.Hsapiens.v75)
library(tidyverse)
library(hdf5r)
library(biovizBase)

パッケージの読み込みをします。

SignacはSeuratと連携して動作し、単細胞アッセイデータ、特にATAC-seqデータの解析を行うためのRパッケージです。 EnsDb.Hsapiens.v75は、Ensemblデータベースのヒトゲノムバージョン75の注釈データを提供するRパッケージです。注釈データには遺伝子、エクソン、イントロン、転写開始部位などのゲノムアノテーション情報を含んでいます。

frag.file <- read.delim('atac_v1_pbmc_10k_fragments.tsv.gz', header = F, nrows = 10)

フラグメントファイルの中身を確認するためのコードです。以下のような出力を得ることができます。

5つのカラムはそれぞれ

V1: 染色体
V2: 染色体上の開始位置
V3: 染色体上の終了位置
V4: バーコード配列
V5: リードペアの総数を表しています。

より詳しい解説は（https://support.10xgenomics.com/single-cell-atac/software/pipelines/latest/output/fragments）をご参照ください。

counts <- Read10X_h5('atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5')

countsの中身は以下のようになっています。

各行は特定の配列を、各列は各細胞を表しています。

chrom_assay <- CreateChromatinAssay(
  counts = counts,
  sep = c(":", "-"),
  fragments = "atac_v1_pbmc_10k_fragments.tsv.gz",
  min.cells = 10,
  min.features = 200
)

SeuratのCreateChromatinAssay関数を用いて、ATAC-seqデータからChromatin Assayオブジェクトを作成するコードです。各オプションについて説明します。

counts = counts:

先ほど作成した疎行列を指定します。

sep = c(“:”, “-“): カウントデータの行名を分割する際の区切り文字を指定します。例: 行名が chr1:713460-714823 のような形式の場合、: と – で分割され、それぞれ chr1、713460、714823 となります。

fragments = “atac_v1_pbmc_10k_fragments.tsv.gz”: フラグメントファイルを指定します。

min.cells = 10: 最低でも10個の細胞で観測される配列のみを保持します。 min.features = 200:

最低でも200個の配列を持つ細胞を保持します。

metadata <- read.csv(file = 'atac_v1_pbmc_10k_singlecell.csv', header = T, row.names = 1)

metadataの読み込みをします。

metadataの中身は上記のようになっており、総フラグメント数やマッピングされなかったリード数、ミトコンドリアにマッピングされたリード数などの情報を持っています。

このmetadataと先ほど作成したChromatin AssayオブジェクトからSeuratオブジェクトを作成します。

pbmc <- CreateSeuratObject(
  counts = chrom_assay,
  meta.data = metadata,
  assay = 'ATAC'
)

CreateSeuratObject関数を使用して、Chromatin AssayオブジェクトとmetadataからSeuratオブジェクトを作成します。それぞれのオプションを説明します。

counts = chrom_assay:
各ピーク（ゲノム領域）と細胞に対するリードカウントが含まれているChromatin Assayオブジェクトを指定します。

meta.data = metadata:
メタデータを指定します。

assay = ‘ATAC’:
解析に使用するアッセイタイプを指定します。今回は、’ATAC’が指定されており、これは単一細胞ATAC-seqデータであることを示しています。

pbmc@assays$ATAC@annotation
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v75)
seqlevels(annotations) <- paste0('chr', seqlevels(annotations))

“annotation”の出力がこちらです。最初は染色体を示す”Rle”列はX,MT等の書式で書かれています。これをseqlevels関数を用いて”chrX”, “chrMT”のような形式に変更します。コード実行後は”annotation”の”Rle”列が変更されているのがわかります。

Annotation(pbmc) <- annotations
pbmc@assays$ATAC@annotation

先ほど作成した”annotation”を”pbmc”に紐づけします。

ここまででSeuratオブジェクトに必要なデータの統合が完了しました。以下からQC作業に移ります。

pbmc <- NucleosomeSignal(pbmc)

NucleosomeSignal関数を用いてヌクレオソームシグナルの計算を行います。

ヌクレオソームシグナルとは

ヌクレオソームに含まれる断片/ヌクレオソームに含まれない断片

で計算されます。一般的にヌクレオソームシグナルが大きい→ヌクレオソームに含まれる断片が多く、ヌクレオソームフリーな断片が少ないことを示します。反対にヌクレオソームシグナルが大きいことはヌクレオソームフリー領域の断片が多いことを示します。

pbmc <- TSSEnrichment(object = pbmc, fast = FALSE)

この関数は、ATAC-seqデータを用いて各細胞の転写開始部位（TSS）周辺のクロマチンアクセスビリティを計算します。特に、TSS周辺のクロマチン状態を評価することで、クロマチンアクセスビリティの全体的な品質を判断できます。

object = pbmc: Seuratオブジェクト pbmc を指定しています。 fast = FALSE: 計算方法の選択肢で、FALSE は通常の計算速度での処理を指定します。高速計算を行いたい場合は TRUE を使用できますが、FALSE の方が精度が高い計算が行われます。

pbmc$blacklist_ratio <- pbmc$blacklist_region_fragments / pbmc$peak_region_fragments
pbmc$pct_reads_in_peaks <- pbmc$peak_region_fragments / pbmc$passed_filters * 100

1行目のコードは、ブラックリスト領域にマッピングされたフラグメントの比率を計算し、新しいメタデータ項目 blacklist_ratio として保存します。ブラックリスト領域とは、ゲノムの中で解析において信頼性が低い領域を指します。

pbmc$blacklist_region_fragments: ブラックリスト領域にマッピングされたフラグメントの数。

pbmc$peak_region_fragments: ピーク領域（クロマチンアクセスビリティが高いとされる領域）にマッピングされたフラグメントの数。

これらを利用してblacklist_ratioを計算します。

2行目のコードは、フィルタリングを通過したフラグメントのうち、ピーク領域にマッピングされたフラグメントの割合を計算し、新しいメタデータ項目 pct_reads_in_peaks として保存します。

pbmc$peak_region_fragments: ピーク領域にマッピングされたフラグメントの数。

pbmc$passed_filters: 全体のフラグメントのうち、フィルタリングを通過したフラグメントの総数。

これらを利用してpct_reads_in_peaks を計算します。

a1 <- DensityScatter(pbmc, x = 'nCount_ATAC', y = 'TSS.enrichment', log_x = TRUE, quantiles = TRUE)
a2 <- DensityScatter(pbmc, x = 'nucleosome_signal', y = 'TSS.enrichment', log_x = TRUE, quantiles = TRUE)

DensityScatter関数を使用して密度散布図を作成します。出力として以下の図が得られます。

左図は横軸にATAC-Seqリード数、縦軸にTSS enrichmentをプロットしています。

右図は横軸にヌクレオソームシグナルを、縦軸にTSS enrichmentをプロットしています。

図中の赤い4本の線はそれぞれ5,10,90,95%点を表しています。

VlnPlot(object = pbmc, 
        features = c('nCount_ATAC', 'nFeature_ATAC', 'TSS.enrichment', 'nucleosome_signal', 'blacklist_ratio', 'pct_reads_in_peaks'),
        pt.size = 0.1,
        ncol = 6)

バイオリンプロットを描画します。以下が出力です。

pbmc <- subset(x = pbmc,
               subset = nCount_ATAC > 3000 &
                 nCount_ATAC < 30000 &
                 pct_reads_in_peaks > 15 & 
                 blacklist_ratio < 0.05 &
                 nucleosome_signal < 4 &
                 TSS.enrichment > 3)

これまでの出力結果を基にフィルタリングを行うコードです。

具体的な内容は以下になります。

nCount_ATAC > 3000: 各細胞におけるリードカウントの合計数が3000以上の細胞を残す。
nCount_ATAC < 30000: リードカウントが30000未満の細胞を残す。
pct_reads_in_peaks > 15: 全リードのうち、ピーク領域内に存在するリードの割合が15%以上の細胞を残す。
blacklist_ratio < 0.05: ブラックリスト領域にマッピングされたリードの割合が5%未満の細胞を残す。
nucleosome_signal < 4: ヌクレオソームシグナルが4未満の細胞を残す。
TSS.enrichment > 3: TSS（転写開始部位）周辺のエンリッチメントスコアが3以上の細胞を残す。

pbmc <- RunTFIDF(pbmc) # normalization
pbmc <- FindTopFeatures(pbmc, min.cutoff = 'q0') # selecting top features
pbmc <- RunSVD(pbmc) # dimensionality reduction

DepthCor(pbmc)

RunTFIDF関数を用いて、TF-IDF（Term Frequency-Inverse Document Frequency）正規化を行い、各ピークの相対的な重要性を計算します。この正規化方法は、各ピークの出現頻度を他のピークの出現頻度と比較し、珍しいピークに対して高いスコアを付与します。

次にFindTopFeatures関数を用いて、最も重要な特徴量を選択します。

min.cutoff = ‘q0’ は、全てのピークを含める設定です。具体的な数値や四分位数（例えば q75 など）を指定することで、上位の特徴量のみを選択することも可能です。

RunSVD関数を用いて、特異値分解による次元削減を行い、データの主成分を抽出します。

DepthCor関数を用いて、データ深度と各成分との相関を計算します。ここでのデータ深度とは、各細胞における総リード数を指します。

このコードを実行すると以下の図が得られます。

図の縦軸はデータ深度（各細胞の総リード）と各成分の値の相関係数を表しており、横軸は次元削減成分を表しています。

グラフを見る限り、多くの成分とデータ深度との間には中程度から強い負の相関が見られます。特に、成分1では顕著な負の相関が見られます。

強い負の相関が見られる成分は、データの深さによる影響を受けている可能性があります。つまりリード数が多いほどその成分の値が小さくなる傾向があるため、バッチ効果や技術的な変動によりデータの均質性が損なわれている可能性があります。そのため、以降は成分1は除いたうえで解析を行います。

pbmc <- RunUMAP(object = pbmc, reduction = 'lsi', dims = 2:30)
pbmc <- FindNeighbors(object = pbmc, reduction = 'lsi', dims = 2:30)
pbmc <- FindClusters(object = pbmc, algorithm = 3)

DimPlot(object = pbmc, label = TRUE) + NoLegend()

UMAP（Uniform Manifold Approximation and Projection）をRunUMAP関数で実行し、データをさらに低次元空間に埋め込みます。これにより、データのクラスタリングパターンをより明確に視覚化できます。最後にDimPlot関数を使って、細胞のクラスタリングを視覚化します。結果として得られたグラフを以下に示します。