S/N比の評価手法 その2 Cross-correlation profile
これは前回の記事の続きですので、未読の方はそちらから読んでください。
Cross-correlation profile
Cross-correlation profile (以下CCP)はENCODEのグループによって提案されたS/N比計測手法です*1。 CCPは、順鎖ー逆鎖間のマップリード分布の「ずれ」を利用します。
この図で示すように、シーケンサはChIPで得られたDNA断片の片端を固定長で読むため、順鎖と逆鎖にマップされるリード分布にはDNA断片長（約150塩基対）ぶんだけ「ずれ」が発生します。得られるピークも順鎖と逆鎖でずれが起きますので、全てのピーク抽出ツールはこのずれを内部で補正しています。
ということは、この「ずれ」の長さだけ順鎖と逆鎖をずらしてやった時にリード分布が重なるはずです。やってみましょう。順鎖と逆鎖をずらしながら、互いのリード分布の相関を取ったものが下図です。
これがCCPです。横軸は順鎖と逆鎖をずらした長さ、縦軸がリード分布の相関です。期待通り、DNA断片長（赤点線）だけずらしたところで相関が最大になっていることがわかります。通常、DNA断片長はpaired-endで読まなければ計測できませんが、CCPを用いることでsingle-end データからでもDNA断片長を推定することが可能です。
青点線はリード長を表します。CCPを描画するとリード長部分にもスパイクが現れることが経験的に知られていますが、これはリピートにマップされるリードの量が影響しています*2。
CCPを作成
Phantompeakqualtoolsのインストール
では、実際にCCPを描いてみましょう。CCPはPhantompeakqualtoolsというツールを用いて計測することができます。 本家サイト とGitHubがありますが、GitHubの方がバージョンが新しいので、そちらを使います。（本家の方はインストールでエラーが出るかも）
git clone https://github.com/kundajelab/phantompeakqualtools.git
cd phantompeakqualtools
R CMD INSTALL spp_1.14.tar.gz
phantompeakqualtools は内部でsppというピーク抽出ツールを利用しています。sppはMACSと並んで最も初期から使われているピーク抽出プログラムのひとつで、R上で動きます。
phantompeakqualtoolsディレクトリの中にあるsppのソースコード (spp_1.14.tar.gz) がありますので、R CMD INSTALLを使ってインストールします。インストールに管理者権限が必要な場合は以下のようにsudoをつけてください。
sudo R CMD INSTALL spp_1.14.tar.gz
データダウンロード
続いてChIP-seqデータをダウンロードします。データはBAMファイルであれば何でもいいのですが、ここではHeLaS3のCTCFとInputを使います。
wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeUwTfbs/wgEncodeUwTfbsHelas3CtcfStdAlnRep1.bam wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeUwTfbs/wgEncodeUwTfbsHelas3InputStdAlnRep1.bam
Phantompeakqualtoolsの実行
ダウンロードが完了したら、phantompeakqualtoolsディレクトリにある run_spp.R を実行しましょう。
# BAMファイル名を変数bamに格納 bam=wgEncodeUwTfbsHelas3CtcfStdAlnRep1.bam # run_spp.R を実行 Rscript run_spp.R \ -p=4 \ -c=$bam -i=$bam \ -savn=Helas3Ctcf.narrowPeak \ -savr=Helas3Ctcf.regionPeak \ -out=Helas3Ctcf.resultfile \ -savp=Helas3Ctcf.pdf # 以下はInputサンプル用 bam=wgEncodeUwTfbsHelas3InputStdAlnRep1.bam Rscript run_spp.R \ -p=4 \ -c=$bam -i=$bam \ -savn=Helas3Input.narrowPeak \ -savr=Helas3Input.regionPeak \ -out=Helas3Input.resultfile \ -savp=Helas3Input.pdf
sppはピーク抽出ツールのため、コマンドもピーク抽出っぽい感じになっています。
-cにChIPサンプル、-iにInputサンプルを指定しますが、やりたいことはS/N比の計測なので、ここでは-cと-i両方に同じBAMファイルを指定しています（つまり対応するInputサンプルは必要ありません）。
-p は使うCPU数ですので、適宜変更してください。Rのマルチスレッディングの環境によっては、-pを指定すると途中でハングして完了しないことがあるので、その場合は-pオプションを削除してシングルCPUで実行してください。
"narrowPeak", "regionPeak"はピークリストで、"resultfile"にS/N比のスコアが記載されます。”.pdf"ファイルが描画されたCCPです。
CCPの実行には少し時間がかかります。-cと-iに同じサンプルを入れているため抽出されるピーク数は当然0です。そのため以下のような警告が出ますが、気にしなくてOKです。
calculating statistical thresholds FDR 0.01 threshold= Inf Detected 0 peaks 警告メッセージ: min(npld$y[npld$fdr <= fdr]) で: min の引数に有限な値がありません: Inf を返します
CCP表示
では、CTCFサンプルのCCPを見てみましょう。
ちゃんとDNA断片長のところが最大になっています。リード長部分のスパイクはほとんど見えないので、これはよいサンプルですね。
下に表示されている"strand-shift (130)"が推定されたDNA断片長です。約130塩基対ということなので、やや短めです。
NSC, RSCがS/N比のスコアで、NSCはバックグラウンドに対するDNA断片長（赤点線）部のピーク高さ、RSCはリード長（青点線）部に対するDNA断片長（赤点線）部のピーク高さを表します。
QtagはNSCとRSCの値を元に判定されるサンプルクオリティの（S/N比）のレベルを表しています。-2,-1,0,1,2の５つの値を持ち、2が最高、-2が最低です。このサンプルは2なので、最高レベルです。
続いてInputサンプルの方を見てみましょう。
Inputサンプルはピークがほとんど存在しないので、DNA断片長部分にスパイクはありません。一方リード長部分に大きなスパイクがあるので、リピート由来のリードをそれなりに含んでいることが予想されます。Qtagは-1ですが、Inputサンプルの場合は当然S/N比が低い方がよいので、Qtagは低い方が望ましいです。
CCPの長所・短所
CCPの長所は、得られるスコアがピーク抽出に依存していないため、サンプルのリード数やピーク抽出法の影響を気にせず利用できる点です。このアドバンテージは大きく、ENCODEやROADMAPではCCPによる品質評価が公式に採用されています。
一方CCPはSharpなピークを持つサンプルにしか使えないという弱点があります。BroadなピークやFaire-seqサンプルなどはピークの強度が弱いためDNA断片長部分に強いスパイクが出にくく、結果として「S/N比の低い悪いサンプル」と判定されてしまいます。そのためCCPの開発者も、「このスコアが悪いことがただちに解析に使えないことを意味するのではない」と注釈を入れています。 私が参加していたIHECプロジェクトでは当初CCPを採用していましたが、後に「Broadサンプルも含めて評価できる」として提案されたdeepToolsに変更されました。
という訳で、次回deepTools、その次にSSPを解説しようと思います。長くなりますね…。
