Model-based Analysis of ChIP-Seq (MACS)

• 著者: Yong Zhang, Tao Liu, Clifford A. Meyer, Jérôme Eeckhoute, David S. Johnson, Bradley E. Bernstein, Chad Nusbaum, Richard M. Myers, Myles Brown, Wei Li, X. Shirley Liu
• Corresponding author: Wei Li (Baylor College of Medicine); X. Shirley Liu (Dana-Farber Cancer Institute / Harvard School of Public Health)
• 雑誌: Genome Biology
• 発行年: 2008
• Epub日: 2008-09-17
• Article種別: Original Article
• PMID: 18798982

背景

Chromatin Immunoprecipitation (免疫沈降) 技術と次世代シーケンシング技術を組み合わせた ChIP-seq (chromatin immunoprecipitation sequencing) は、ゲノムワイドなタンパク質-DNA相互作用を同定する強力な手法である cistrome (ゲノムワイドな転写因子結合領域の総体) 解析として、従来の ChIP-chip (chromatin immunoprecipitation on chip) 技術を急速に置換しつつあった。先行研究では、ゲノムタイルマイクロアレイを用いた ChIP-chip がゲノムワイドな転写因子結合部位の同定に貢献してきたが、解像度やコストの面で限界があった。その後、高スループットなシーケンシング技術の登場により、より高精度なマッピングが可能となった。しかし、初期の ChIP-seq 解析においては、いくつかの重大な技術的課題が残されていた。第一に、シーケンスされるリード (tag) は免疫沈降された DNA 断片の両端に過ぎず、実際のタンパク質結合部位から数十から数百塩基対 (bp) 離れた位置にマッピングされる。この tag から結合部位へのシフトサイズは実験ごとに異なり、ユーザーにとっては未知のパラメータであった。第二に、ゲノム上の tag 分布には、局所的なクロマチン構造、シーケンスバイアス、マッピングバイアス、およびゲノムコピー数変動に起因する顕著な偏りが存在することが明らかになっていた。当時公表されていたいくつかの ChIP-seq 研究のうち、Barski et al. (2007) Barski et al などのデータセットでは、コントロールサンプルが十分に活用されておらず、局所的な背景ノイズを補正する標準的な統計モデルが確立されていなかった。特に、シグナルが弱い領域ではコントロールの tag 数が極めて少なく、局所的な背景バイアスを頑健に推定するためのアプローチが不足していた。このように、ChIP-seq データの特性に最適化された、バイアス補正と高解像度なピーク検出を両立するアルゴリズムは未確立であり、解析手法の標準化に向けた大きな gap が残されていた。さらに、Johnson et al. (2007) Johnson et al や Lupien et al. (2008) Lupien et al などの研究でも、シーケンスデータの特性を最大限に活かすための数理モデルが不足しており、実験バイアスを排除しつつ真の結合シグナルを正確に検出する手法の開発が強く望まれていた。このように、局所的なゲノムバイアスを動的に補正する手法は未確立であり、解析の信頼性を担保する統計モデルが不足していた。

目的

本研究の目的は、ChIP-seq データ解析における上記の問題点を克服し、高精度かつ高分解能なピーク検出を実現する新しいアルゴリズム「MACS (Model-based Analysis of ChIP-Seq)」を開発することである。具体的には、以下の3つの技術的目標を達成することを目指した。第一に、Watson 鎖と Crick 鎖にマッピングされた tag の双峰性 (bimodal) パターンを利用して、結合部位へのシフトサイズ (d) を経験的かつ自動的にモデル化し、予測される結合位置の空間分解能を大幅に向上させること。第二に、コントロールサンプルの局所的な tag 密度を動的に評価する「動的ポアソン分布 (dynamic Poisson distribution)」を導入し、ゲノム全体のバイアスや局所的なノイズを効果的に補正すること。第三に、既存のピーク検出アルゴリズム (FindPeaks, ChIPSeq Peak Finder, QuEST) と比較して、感度、特異性、空間分解能のすべての面において MACS の優位性を実証し、コントロールサンプルの有無に関わらず頑健な解析を可能にすることである。これにより、次世代シーケンシング技術を用いたゲノムワイドなタンパク質-DNA相互作用解析の標準的なプラットフォームを確立することを目指した。

結果

**双峰性パターンに基づくシフトサイズ d の正確な推定**: MACS は、ChIP-seq データにおける Watson 鎖と Crick 鎖の tag が示す双峰性 (bimodal) パターンを利用して、結合部位へのシフトサイズ d を経験的にモデル化する。MCF-7 細胞株を用いた FOXA1 ChIP-seq データ (モデル構築のために n=1000 model peaks をサンプリング) に適用したところ、ソニケーションサイズが約 500 bp、サイズ選択が約 200 bp の条件下において、MACS は d=126 bp (すなわち 3’ へのシフトサイズは 63 bp) と推定した (Fig 1)。この推定値は、実際の FKHR (forkhead in rhabdomyosarcoma) モチーフ配列のアライメントから直接計算された d=122 bp と極めて近い値であった。同様に、NRSF データセットでは MACS による推定値 d=96 bp に対し、NRSE (neuron-restrictive silencer element) モチーフに基づく推定値は 70 bp であり、CTCF データセットでは MACS による推定値 d=76 bp に対しモチーフに基づく推定値は 62 bp であった。このように、MACS は事前知識なしに、tag のゲノム分布のみから高精度にタンパク質-DNA結合部位の位置を補正できることが実証された。

**動的ポアソン分布による局所バイアスの排除と偽陽性の抑制**: ゲノム上の局所的なシーケンスバイアスやコピー数変動を補正するため、MACS は動的パラメータ λlocal を導入している。FOXA1 の ChIP サンプルとコントロールサンプルの比較において、両者の tag 数はゲノム全体で強く相関していた (Fig 1)。一様なバックグラウンド分布 λBG のみを用いる従来のモデルでは、コントロールサンプルでも tag 数が局所的に増加している領域において、多くの偽陽性ピークが検出されてしまう。MACS の λlocal は、これらのバイアス領域において λBG の 2 から 10 倍に動的に上昇し、偽陽性を効果的に排除した。コントロールサンプルが存在する場合、n=7000 peaks を検出した際の偽発見率 (FDR) は 0.4% と極めて低く抑えられた。さらに、コントロールサンプルが存在しない場合を想定し、ChIP サンプル自身から λlocal を推定した場合でも、FDR は 3.8% に維持された。これに対し、一様な λBG を用いた場合の FDR は 41.2% まで悪化しており、動的ポアソン分布モデルがコントロールなしの解析においても極めて頑健に機能することが示された。

**既存アルゴリズムとの比較における感度と空間分解能の優位性**: MACS の性能を、既存のピーク検出ツールである ChIPSeq Peak Finder、FindPeaks、および QuEST と比較評価した。FOXA1 データセットにおいて、MACS は p<10⁻⁵ のカットオフで n=3449 peaks を検出し、そのうち 87% のピークが中心から 50 bp 以内に既知の FKHR モチーフを含んでいた (Fig 2)。これは既存ツールと同等以上のモチーフ含有率でありながら、より多くの真のピークを検出していることを意味する。また、ピークの頂点 (summit) から最も近いモチーフ中心までの平均距離を測定したところ、MACS は中央値 11 bp という最も短い距離を達成し、優れた空間分解能を示した。NRSF データセットにおいても、MACS は n=1945 peaks を検出し、その 72% に NRSE モチーフが含まれており、既存手法よりも少ない偽陽性数で高い感度を達成した (Fig 2)。

**ChIP-Seq と ChIP-chip プラットフォーム間の検出性能比較**: FOXA1 を対象として、ChIP-seq (MACS 検出、FDR <1%) と従来の ChIP-chip (MAT (Model-based Analysis of Tiling-arrays) 検出、FDR <1%, fold change >2.0) の比較を行った。その結果、ChIP-seq で検出されたピークの 65.4% が ChIP-chip でも検出されており、両プラットフォーム間で高い整合性が確認された (Fig 3)。ChIP-seq のみで検出されたユニークなピーク (21.4%) のうち、n=1684 peaks について ChIP-chip のシグナル強度を再確認したところ、閾値以下ではあるものの、ゲノムバックグラウンドと比較して有意に高いシグナルが観察された (Wilcoxon 検定、p<10⁻³²⁰)。また、ChIP-seq で得られたピークの幅は ChIP-chip の約半分であり、モチーフ中心への近接度も有意に優れていた (Fig 3)。これにより、ChIP-seq が従来のマイクロアレイ技術と比較して、より高い空間分解能と検出感度を有していることが統計的に証明された。

**シーケンス深度と飽和度に関する検証**: MACS は、ユーザーがシーケンス深度の十分性を評価できるように、異なるサブサンプリング比率 (90% から 20% の tag を使用) におけるピーク検出数をシミュレーションする飽和度評価機能を備えている。FOXA1 (3.9 million tags) および NRSF (2.2 million tags) のデータを用いた検証により、60-fold 以上の高い濃縮度を持つピークは比較的少ないリード数で飽和に達する一方、40-fold 以下の弱いシグナル領域を完全に網羅するためには、さらなるシーケンス深度が必要であることが明らかになった (Fig 4)。この機能により、実験デザインに応じた最適なシーケンス深度の決定が可能となった。

考察/結論

先行研究との違い: MACS は、従来の ChIP-chip 解析手法や初期の ChIP-seq ピーク検出アルゴリズムと異なり、tag の双峰性パターンからシフトサイズ d を動的に推定し、すべての tag 位置を 3’ 方向に d/2 だけシフトさせるモデルベースのアプローチを採用している。従来の多くの手法が固定されたウィンドウサイズや単純なスライディングウィンドウによるカウントに依存していたのに対し、MACS はシーケンスデータの特性自体から物理的な断片サイズを経験的に逆算する。また、一様なバックグラウンド分布を仮定する先行手法とは対照的に、局所的なゲノムバイアスを 1 kb、5 kb、10 kb のマルチスケールウィンドウで評価する動的ポアソン分布モデルを導入した点が決定的に異なる。

新規性: 本研究は、ChIP-seq データに特有の局所的バイアス (クロマチン構造やシーケンスバイアス) を、コントロールサンプルの有無に関わらず頑健に補正できるアルゴリズムを新規に開発した。特に、コントロールサンプルが存在しない場合でも、ChIP サンプル自体の局所情報から動的パラメータ λlocal を推定し、FDR を 3.8% という極めて低い水準に抑える手法は、本研究で初めて提案されたものである。これにより、高価なシーケンスコストを抑えつつ、信頼性の高い cistrome 解析を行うための新しい標準的枠組みが確立された。

臨床応用: 本アルゴリズムの臨床的意義は極めて大きい。がんゲノム研究やエピゲノミクス解析において、転写因子やヒストン修飾の結合プロファイルを正確に把握することは、治療標的の同定やバイオマーカー開発に直結する。例えば、肺がんにおける FOXA1 や NKX2-1 (TTF-1) の cistrome 解析、あるいは小細胞肺がん (SCLC; small cell lung cancer) のサブタイプ分類 (ASCL1, NEUROD1, POU2F3) におけるエンハンサー領域の同定など、臨床現場に直結する基礎研究において MACS は必須のツールとして広く活用されている。疾患特異的なエピゲノム変化をハイスループットかつ高解像度で同定できるため、精密医療 (precision medicine) の推進に貢献している。

残された課題: 一方で、いくつかの limitation や今後の課題も存在する。第一に、MACS は転写因子のような sharp なピーク検出には最適化されているが、H3K27me3 や H3K4me1 などの広範なヒストン修飾領域 (broad domain) の検出には不向きであり、これには後継バージョンである MACS2 (Model-based Analysis of ChIP-Seq version 2) の broad peak モードなどの開発が必要であった。第二に、PCR 重複リードの過度な除去アルゴリズムは、極めて高いシーケンス深度において真のシグナルクラスターを誤って排除するリスクがある。第三に、ChIP-exo や CUT&RUN、CUT&Tag といった超高解像度技術の登場に伴い、塩基レベルの解像度 (base-pair resolution) を求める解析においては、MACS の空間分解能には限界がある。これらの課題に対する今後の研究とアルゴリズムの改良が期待される。

方法

MACS アルゴリズムの実装および評価は、以下の手順に従って行われた。 アルゴリズムは大きく4つのステップで構成される。 Step 1 (シフトサイズ d のモデル化): ユーザーが指定したソニケーションサイズ (bandwidth) と、高信頼な濃縮度 (mfold) をパラメータとして用いる。ゲノム全体にわたって 2×bandwidth のウィンドウをスライドさせ、ランダム分布から mfold 以上の濃縮を示す領域を探索する。これらの高信頼ピークからランダムに n=1000 model peaks をサンプリングし、Watson 鎖と Crick 鎖の tag 分布をアライメントする。両鎖 of ピークモード間の距離を d と定義し、すべての tag を 3’ 方向に d/2 だけシフトさせることで、実際のタンパク質-DNA結合部位に位置を補正する。 Step 2 (動的ポアソン分布を用いたピーク検出): ゲノム全体のバックグラウンド tag 数から平均パラメータ λBG を計算する。各候補ピークにおいて、局所的なバイアスを考慮するため、1 kb、5 kb、10 kb のウィンドウサイズにおける局所平均 λ1k, λ5k, λ10k を算出する。動的パラメータ λlocal を λlocal = max(λBG, λ1k, λ5k, λ10k) として定義し、Poisson 分布に基づいて各候補ピークの p-value を計算する。デフォルトの閾値は p<10⁻⁵ とした。 Step 3 (偽発見率 FDR (false discovery rate) の推定): ChIP サンプルとコントロールサンプルのラベルを入れ替える「サンプルスワップ (sample swap)」を行い、同様のパラメータで偽ピーク数を算出する。これにより、各 p-value 閾値における経験的 FDR を算出する。 Step 4 (重複 tag の除去): PCR 増幅バイアス等によるアーティファクトを排除するため、同一ゲノム座標にマッピングされた重複 tag 数を制限する。制限値は、シーケンス深度に基づく binomial 分布 (p<10⁻⁵) を用いて動的に決定される。

アルゴリズムの評価には、以下の n=3 cell lines から得られたデータセットを用いた。 (a) MCF-7 乳がん細胞株 (MCF-7 cell line) を用いて本研究で新規に生成した FOXA1 (forkhead box A1) ChIP-seq データ (3.9 million のユニークマッピングリード、および 5.2 million のコントロールリード)。FOXA1 はエストロゲン受容体の結合を先導するパイオニア因子として知られており Lupien et al、ライブラリ調製には、3つの独立した実験から得られた DNA の等量混合物 (n=3 independent experiments) を用いた。 (b) Jurkat T 細胞株 (Jurkat cell line) を用いた NRSF (neuron-restrictive silencer factor) ChIP-seq データ (2.2 million の ChIP リード、2.8 million のコントロールリード) Johnson et al。 (c) CD4⁺ T 細胞 (CD4+ T cells) を用いた CTCF (CCCTC-binding factor) ChIP-seq データ (2.9 million の ChIP リード) Barski et al。なお、CTCF は CCCTC (cytosine-cytosine-cytosine-thymidine-cytosine) 結合因子を指す。 統計的比較には Wilcoxon 符号付き順位検定 (Wilcoxon signed-rank test) を用いて、ピーク幅や空間分解能の有意差を評価した。