Palmsonntagmorgen

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

S/N比の評価手法 その4 SSP

時間がかかってしまいましたが、やっとSSPの登場です。

この記事は以下の記事の続きです。

S/N比の評価手法 その1 - Palmsonntagmorgen

S/N比の評価手法 その2 Cross-correlation profile - Palmsonntagmorgen

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

SSPとは

SSP (Strand-shift profile)は Cross-correlation profile (上記その2参照)を拡張したもので、相関係数を取るかわりに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の描画

その2 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 f:id:rnakato:20180427132049p:plain

Input.jaccard.pdf f:id:rnakato:20180427132126p:plain

左図はCCPとほぼ同じスケールですが、CCPと同じようなプロファイルが得られていることがわかります。 右図はx軸(順鎖と逆鎖のずれ)のスケールを1Mbpまで伸ばした時のプロファイル(log scale)です。
(CCPでは計算量の問題でx軸が1kbp程度、5bpおきの計測となっていますが、SSPは計算が高速なためsingle bp resolutionで遠くまで計算することができます)

左下にスコアが出ていますが、NSCS/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 

タブ区切りの表になっているのですが、列の表示がずれていたり、有効桁数表示が多かったりして見づらいので(やたら桁数表示が多いツールってセンスがないですよね!)、エクセルに貼りつけて棒グラフにしてみましょう。

ぱっ!

f:id:rnakato:20180427145753p:plain:w400

f:id:rnakato:20180427145807p:plain:w400

f:id:rnakato:20180427145826p:plain:w400

ポイントをまとめると、

  • SSPはH3K9me3, H3K4me3ともInputより高いスコアを出しています。リード数にも依存していません。
  • これに対してCCPはH3K9me3とInputの差がほとんどありません。これが「sharp peakにしか使えない」と言われる所以です。
  • deepToolsはSSPと同じくH3K9me3, H3K4me3がInputより高いスコアを出していますが、リード数に依存してスコアが変わってしまっており、スコアがリード数に依存していることがわかります。
    • (2つしかないとわかりにくいですが、1000万リードよりも少なくなると急激に値が減少します)。

まとめ

以上、簡単ですが、SSPが他のツールよりも優れている点を比較評価してみました。 今回は3サンプルのみの簡単な実験ですが、論文では1000以上のChIP-seqサンプルを使った評価をしているので、興味があれば参照してください。

Sensitive and robust assessment of ChIP-seq read distribution using a strand-shift profile | Bioinformatics | Oxford Academic