S/N比の評価手法 その4 SSP
時間がかかってしまいましたが、やっとSSPの登場です。
この記事は以下の記事の続きです。
S/N比の評価手法 その１ - Palmsonntagmorgen
S/N比の評価手法 その2 Cross-correlation profile - Palmsonntagmorgen
S/N比の評価手法 その3 deepTools - Palmsonntagmorgen
SSPとは
SSP (Strand-shift profile)は Cross-correlation profile (上記その２参照)を拡張したもので、相関係数を取るかわりにJaccard indexを使って順鎖と逆鎖のマップリード分布の一致度を計算します。
Jaccard indexを用いるメリットとして以下が挙げられます。
- 相関係数に比べて計算が高速。
- PCR biasをフィルタした状態ではマップリード数の値が0か1しかないため、相関係数は不適。
- unmap-unmapとmap-mapを両方評価する相関係数に対し、Jaccard indexはmap-mapペアのみを評価する。
SSPもいくつか機能があるのですが、今回はS/N比の計測機能のみに注目します。
SSPインストール
aptで必要なライブラリ群をインストールした後、ソースコードをGitHubからダウンロードの後コンパイルするとsspの実行ファイルが作成されます。
$ sudo apt install git build-essential libboost-all-dev \ libgsl-dev libz-dev samtools $ git clone https://github.com/rnakato/SSP.git $ cd SSP $ make
以下のコマンドで実行ファイルにパスを通します。
(SSPをダウンロードしたディレクトリが $HOME/my_chipseq_exp/SSP の場合)
参考:
環境変数PATHの通し方 - Palmsonntagmorgen
$ export PATH = $PATH:$HOME/my_chipseq_exp/SSP/bin
sspとタイプしてヘルプが表示されればインストール成功です。
$ ssp 11:30:18 SSP v1.1.1 =============== Usage: ssp [option] -i <inputfile> -o <output> --gt <genome_table> Use --help option for more information on the other options
Strand-shift profileの描画
その２ CCP の時と同じCTCFとInputのサンプルを使って、Strand-shift profileを描画してみましょう。
$ bam=wgEncodeUwTfbsHelas3CtcfStdAlnRep1.bam $ ssp -i $bam -o CTCF -p 4 --gt genometable.hg19.txt --mptable mptable.hg19.txt $ bam=wgEncodeUwTfbsHelas3InputStdAlnRep1.bam $ ssp -i $bam -o Input -p 4 --gt genometable.hg19.txt --mptable mptable.hg19.txt
-i で入力サンプルを指定します。BAM, SAM, tagAlign, Bowtieフォーマットが入力可能で、ファイルのフォーマットはファイル名から自動で判定されます。
--gt はgenome tableファイル、 --mp はマップ可能な(mappable) ゲノム領域の長さを記したファイルです。いずれもインストールされたSSPフォルダのdata ディレクトリの中にあります。-pは使用するCPU数です。
出力ファイルはsspoutディレクトリの下に生成されます。CTCF.jaccard.pdf, Input.jaccard.pdf を見ると以下のようになっています。
CTCF.jaccard.pdf
Input.jaccard.pdf
左図はCCPとほぼ同じスケールですが、CCPと同じようなプロファイルが得られていることがわかります。
右図はｘ軸（順鎖と逆鎖のずれ）のスケールを1Mbpまで伸ばした時のプロファイル（log scale）です。
(CCPでは計算量の問題でx軸が1kbp程度、5bpおきの計測となっていますが、SSPは計算が高速なためsingle bp resolutionで遠くまで計算することができます)
左下にスコアが出ていますが、NSCがS/N比の値として用いられるスコアです。CTCFは42.9, Inputは1.42 となっており, S/N比の差を正しく反映していそうです。
SSP, CCP, deepToolsの比較評価
それでは、SSP, CCP, deepToolsから得られるS/N比のスコアを比較してみましょう。
ここでの実験では以下のようにします。
- sharp peak, broad peak, input の三種類のサンプルを用いる。
- 1000万リード、2000万リードをランダムに抽出したファイルをそれぞれ作り、スコアが同じになるか比較（read depthの影響）。
サンプル数が多くなりますので、シェル変数やfor文を利用していきます。シェルスクリプトに不慣れな方は以下の記事も参考にしてください。
DROMPA3: その5 シェル変数を使う - Palmsonntagmorgen
データダウンロード
データはROADMAP projectのE018（iPS-15b cells）から、H3K4me3 (sharp), H3K9me3 (broad), Inputをダウンロードします。
for prefix in E018-H3K4me3 E018-H3K9me3 E018-Input; do wget http://egg2.wustl.edu/roadmap/data/byFileType/alignments/consolidated/$prefix.tagAlign.gz done
データ変換
ダウンロードされたサンプルは3000万マップリードありますが、ここから1000万、2000万リードを抽出したファイルを生成します。生成には SSP/scripts フォルダにある randompickup4bed.pl を使います。
for prefix in E018-H3K4me3 E018-H3K9me3 E018-Input; do file=$prefix.tagAlign.gz gunzip $file # gzファイルを解凍 $prefix.tagAlign ができる for num in 10 20; do # $prefix.tagAlign から ${num}000000 リードをランダムに抽出 randompickup4bed.pl $prefix.tagAlign ${num}000000 > $prefix.${num}mil.tagAlign # $prefix.${num}mil.tagAlign を圧縮 $prefix.${num}mil.tagAlign.gzができる gzip $prefix.${num}mil.tagAlign done done
deepToolsは sorted bam 形式しか受け付けませんので、tagAlignファイルをsorted bamに変換します。
for prefix in E018-H3K4me3 E018-H3K9me3 E018-Input; do for num in 10 20; do bedToBam -i $prefix.${num}mil.tagAlign.gz -g genometable.txt | samtools sort -@4 > $prefix.${num}mil.sort.bam samtools index $prefix.${num}mil.sort.bam done done
完了すると、以下のようにファイルができているはずです。
$ ls
E018-H3K4me3.10mil.sort.bam E018-H3K9me3.10mil.sort.bam E018-Input.10mil.sort.bam
E018-H3K4me3.10mil.sort.bam.bai E018-H3K9me3.10mil.sort.bam.bai E018-Input.10mil.sort.bam.bai
E018-H3K4me3.10mil.tagAlign.gz E018-H3K9me3.10mil.tagAlign.gz E018-Input.10mil.tagAlign.gz
E018-H3K4me3.20mil.sort.bam E018-H3K9me3.20mil.sort.bam E018-Input.20mil.sort.bam
E018-H3K4me3.20mil.sort.bam.bai E018-H3K9me3.20mil.sort.bam.bai E018-Input.20mil.sort.bam.bai
E018-H3K4me3.20mil.tagAlign.gz E018-H3K9me3.20mil.tagAlign.gz E018-Input.20mil.tagAlign.gz
E018-H3K4me3.tagAlign E018-H3K9me3.tagAlign E018-Input.tagAlign
実行
生成されたファイルを使ってSSP, CCP, deepToolsを実行します。 かなり時間がかかりますので、気長に待ちましょう。
for prefix in E018-H3K4me3 E018-H3K9me3 E018-Input; do for num in 10 20; do head=$prefix.${num}mil # SSP ssp -i $head.tagAlign.gz -o SSP-$head --gt genometable.txt -p 4 --mptable mptable.txt # CCP Rscript run_spp.R -p=4 -c=$head.tagAlign.gz -i=$head.tagAlign.gz \ -savn=CCP-$head.narrowPeak -savr=CCP-$head.regionPeak \ -out=CCP-$head.resultfile -savp=CCP-$head.pdf # deepTools plotFingerprint -b $head.sort.bam --JSDsample $head.sort.bam --ignoreDuplicates -p 4 \ -plot DT-fingerprint.$head.png --outQualityMetrics DT-$head.txt done done
結果表示
結果を見てみましょう。ひとつひとつファイルを確認するのは時間がかかるので、シェルスクリプトでぱぱっとやってしまいます。 （それぞれの行が何をしているのか考えてみてください）
echo -e "Sample\tSSP\tCCP\tdeepTools" for prefix in E018-H3K4me3 E018-H3K9me3 E018-Input; do for num in 10 20; do head=$prefix.${num}mil echo -en "`tail -n1 sspout/SSP-$head.stats.txt | cut -f1,6`\t" echo -en "`cut -f9 CCP-$head.resultfile`\t" echo -e "`tail -n1 DT-$head.txt| cut -f9`" done done
表示される結果は以下のようになります。
Sample SSP CCP deepTools E018-H3K4me3.10mil 13.1502 1.429187 0.42114196480846744 E018-H3K4me3.20mil 13.332 1.427474 0.4537138959592519 E018-H3K9me3.10mil 3.37982 1.062391 0.2666993695154958 E018-H3K9me3.20mil 3.36639 1.062076 0.29083049547676265 E018-Input.10mil 1.31723 1.012775 0.1975224100365595 E018-Input.20mil 1.29461 1.012816 0.20723845047525988
タブ区切りの表になっているのですが、列の表示がずれていたり、有効桁数表示が多かったりして見づらいので（やたら桁数表示が多いツールってセンスがないですよね！）、エクセルに貼りつけて棒グラフにしてみましょう。
ぱっ！
ポイントをまとめると、
- SSPはH3K9me3, H3K4me3ともInputより高いスコアを出しています。リード数にも依存していません。
- これに対してCCPはH3K9me3とInputの差がほとんどありません。これが「sharp peakにしか使えない」と言われる所以です。
- deepToolsはSSPと同じくH3K9me3, H3K4me3がInputより高いスコアを出していますが、リード数に依存してスコアが変わってしまっており、スコアがリード数に依存していることがわかります。
- （２つしかないとわかりにくいですが、1000万リードよりも少なくなると急激に値が減少します）。
まとめ
以上、簡単ですが、SSPが他のツールよりも優れている点を比較評価してみました。 今回は3サンプルのみの簡単な実験ですが、論文では1000以上のChIP-seqサンプルを使った評価をしているので、興味があれば参照してください。
