Palmsonntagmorgen

NGSデータを使った解析と、その周辺。

S/N比の評価手法　その3 deepTools

ChIP-seq NGS 品質評価

今回はIHEC projectの公式品質評価ツールに採用されているdeepToolsについて解説します。

deepToolsとは

deepToolsはChIP-seq, RNA-seq, MNase-seqなどの品質評価及び種々の可視化をするために作られたソフトウェアで、samtoolsやbamtoolsでは手の届かないようなちょっと込み入った解析をすることができます。 pythonで書かれたツールであり、Anacondaでインストールできるほか、Galaxyブラウザ上で利用することもできます。

色々なお役立ち機能があるのでそのうち個別に利用事例を紹介したいと思いますが、今回はS/N比の計測機能を利用します。

deepToolsインストール

condaを使ってインストールします。deepToolsはPython2系にも3系にもインストール可能です。

conda install -c bioconda deeptools

conda (Anaconda)とは何ぞや？という方は以下のエントリを参考にしてください。

pyenvでPython環境を構築する - Palmsonntagmorgen

データダウンロード

データは前回の記事と同じものを使います。新しくダウンロードする場合は以下のコマンドを実行してください。

wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeUwTfbs/wgEncodeUwTfbsHelas3CtcfStdAlnRep1.bam
wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeUwTfbs/wgEncodeUwTfbsHelas3InputStdAlnRep1.bam

deepToolsはbam indexを要求しますので、ダウンロードしたBAMファイルのindexを作成します。

samtools index wgEncodeUwTfbsHelas3CtcfStdAlnRep1.bam 
samtools index wgEncodeUwTfbsHelas3InputStdAlnRep1.bam

S/N比の計測

S/N比の計測には、plotFingerprintコマンドを使います。

plotFingerprint -b wgEncodeUwTfbsHelas3CtcfStdAlnRep1.bam \
    -plot fingerprint.png --ignoreDuplicates -p 2 \
    --outQualityMetrics CTCF.qc.txt \
    --JSDsample wgEncodeUwTfbsHelas3InputStdAlnRep1.bam

-bでChIPサンプル、--JSDsample でInputサンプルを指定します。plotFingerprintコマンドでのS/N比の計測にはChIPサンプルに対応するInputサンプルが必要になります（詳細は後述）。
-plot fingerprint.png , --outQualityMetrics CTCF.qc.txtでそれぞれ出力される図、スタッツのファイル名を指定します。 --ignoreDuplicates はPCR biasをフィルタするオプションで、通常このオプションをつけることが推奨されています。 -p 2は使用するCPU数です。

fingerprint.plot

得られたfingerprint.png を表示してみましょう。

f:id:rnakato:20180407214948p:plain

この図は、「各binにマップされたリード数の全体に対する割合」を累積値で表示したものです。 x軸はリード数でソートされたbin、y軸が累積リード数の割合を示しています。
仮にリードの分布が完全に一様だとすると、このプロットはまっすぐ対角線上に伸びるはずです。 一方、特定の領域にリード分布が偏っていれば、x軸のrankが右端に行った時点で急激に伸びるはずです。 つまり、S/N比が高いサンプルほど、このプロットは右下側にずれることになります。

上図を見ると、確かにInputサンプル（オレンジ）よりもChIPサンプル（青）の方が右下によっています。 Inputサンプルはx軸が0.2のあたりまでy軸が0ですが、これはリードがマップされていないbinの割合を示しています。ゲノム全体の二割がマップされていません。一方、CTCFサンプルは6割以上のbinにリードがマップされていないようです。これをS/N比が高いと見るか、ゲノムカバー率が低いと見るかは微妙なところですが、この場合はおそらくリード数の不足によるものでしょう。

Jensen-Shannon distance (JSD)

これをスコアで定量的に表したものが"Jensen-Shannon distance (JSD)"です。CTCF.qc.txtをエクセルで開いてみましょう。

f:id:rnakato:20180407220412p:plain

ChIPサンプルとInputサンプルそれぞれについて、上のプロットにもとづいて計算されたスコアが並んでいます。 各スコアの詳細については公式サイトを参照してください。

今回見たいのは "JS Distance" の項です。JSDはカルバック・ライブラー情報量 (Kullback-Leibler divergence)を改良したもので、ChIPサンプルとInputサンプルのリード分布ラインがどのくらい離れているかをの違いを相対エントロピーとして計算したものです。スコアは0～1の値を持ち、今回のCTCFサンプルでは約0.27となっています。

Synthetic JSD

お気づきのようにこの計算にはInputサンプルが必要なので、用いるInputサンプルによって値が変わります。 また、Inputサンプルの "JS Distance" の項はNAになってしまいます。
これに対し、"Synthetic JS Distance"の項ではInputサンプルにもスコアが表示されています。Synthetic JS Distanceは実際のInputサンプルを用いる代わりに、同じリード数で完全一様なInputサンプルを疑似的に生成することでInputなしでJSDを計算できるようにしたものです。バックグラウンドモデルにはPoisson分布が用いられています。

したがって、Synthetic JSDはInputサンプルなしで計算可能です。やってみましょう。以下のように、ChIPとInputの両方にChIPサンプルを指定してみます。
（--JSDsampleオプションを指定しないと出力からJSDの項が消えてしまうため、ここでは同じサンプルを指定しています）

$ plotFingerprint -b wgEncodeUwTfbsHelas3CtcfStdAlnRep1.bam \
$       -plot fingerprint.noInput.png --ignoreDuplicates \
$       -p 2 \
$       --outQualityMetrics CTCF.noInput.qc.txt \
$       --JSDsample wgEncodeUwTfbsHelas3CtcfStdAlnRep1.bam

出力は以下のようになります。

f:id:rnakato:20180407222458p:plain いくつかのスコアはNAになりますが、Synthetic JSDは先ほどと同じ値が得られています。従って、deepToolsでS/N比を計測する時はSynthetic JSDを用いるのがベストでしょう。

JSDの長所・短所

JSDの良い点は、sharp/broadによらずS/N比を計測できるという点です。また、ゲノム全体でなく、一部の領域をランダムに抽出したうえで計算するので、計算時間もそれほどかかりません。 一方、ピークとリピート領域のような濃縮を区別しないので、リピートに大量にマップされているようなサンプルも高いS/N比として検出されます。また、上の例のように、サンプルのゲノムカバー率が低すぎる場合にはうまく計算できません。JSDの最も大きな問題は、スコアがリード数に依存してしまうという点です。たくさんのサンプルをこの手法で比較する場合には、リード数を揃えてやる必要があります。