Palmsonntagmorgen

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

記事一覧

解析環境構築

環境変数PATHの通し方

pyenvでPython環境を構築する

バイオインフォマティクスのためのpython環境構築方法を考える

データ生成

 

常染色体と性染色体のみのゲノム配列ファイル genome.fa を作成する

2bit genome を作成する

genome tableを作成する

 

fastqデータ取得

SRAからfastqを取得する

ENA,DDBJからfastqを取得する

 

ゲノムマッピング

Readをゲノムにマッピング (その1) (2017/12/19 追記あり)

Readをゲノムにマッピング (その2)

Readをゲノムにマッピング (その3) 圧縮ファイルを入力にする方法

 

マップデータ操作

SAMtoolsとリダイレクト

SAMtoolsワンライナー覚書

 

ChIP-seq解析: DROMPA

DROMPA3: その1 インストール

DROMPA3: その2 parse2wig

DROMPA3: その3 ピーク抽出(peak calling)

DROMPA3: その4 マップリード分布の可視化その1

DROMPA3: その5 シェル変数を使う

DROMPA3: その6 ChIP/Input ratio 及び p値の可視化

DROMPA3: その7 -i オプション詳細

DROMPA3: その8 GVコマンドでのマクロな可視化

DROMPA3: その9 リードプロファイル

ChIP-seq解析: 品質評価

Library complexity (PCR bias)とは何か

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

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

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

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

ChIP-seq解析: ピークを入力とする操作

BEDtoolsワンライナー覚書

DROMPA3: その9 リードプロファイル

DROMPA3を用いたリードプロファイルの描画です。
ここでプロファイルと呼んでいるものは、遺伝子回り、あるいはピーク回りにおけるマップリードの平均分布のことで、aggregation plotと呼ばれることもあります。
全遺伝子や全ピークの平均値として見ることにより、通常のピーク抽出では得られない微小なリードの濃縮を捉えたり、発現遺伝子群と非発現遺伝子群の間でのリード濃縮の違いを可視化することができます。

インストール

DROMPA3のインストール方法はこの記事を参照してください。

PROFILEコマンドは内部でRスクリプトを実行しますので、Rが必要になります。

Genome table作成

DROMPA3の実行にはGenome tableファイルが必要になります。 Genome tableの作成は以下の記事を参照してください。

genome tableを作成する - Palmsonntagmorgen

遺伝子アノテーション

可視化する際に必要になる遺伝子アノテーションデータ(refFlat形式)をUCSCのサイトからダウンロードします。

$ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/database/refFlat.txt.gz
$ gunzip refFlat.txt.gz  # 解凍

parse2wig

今回もENCODE projectにあるK562のヒストン修飾群、遺伝子アノテーションを用います。

$ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeUwHistone/wgEncodeUwHistoneK562H3k4me3StdAlnRep1.bam
$ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeUwHistone/wgEncodeUwHistoneK562H3k27me3StdAlnRep1.bam
$ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeUwHistone/wgEncodeUwHistoneK562H3k36me3StdAlnRep1.bam
$ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeUwHistone/wgEncodeUwHistoneK562InputStdAlnRep1.bam
$ parse2wig -i wgEncodeUwHistoneK562H3k4me3StdAlnRep1.bam -o H3K4me3 -gt genometable.txt -f BAM
$ parse2wig -i wgEncodeUwHistoneK562H3k27me3StdAlnRep1.bam -o H3K27me3 -gt genometable.txt -f BAM
$ parse2wig -i wgEncodeUwHistoneK562H3k36me3StdAlnRep1.bam -o H3K36me3 -gt genometable.txt -f BAM
$ parse2wig -i wgEncodeUwHistoneK562InputStdAlnRep1.bam -o Input -gt genometable.txt -f BAM

PROFILEコマンド

drompa_draw PROFILE とタイプするとヘルプが表示されます。

$ drompa_draw PROFILE
Usage: drompa_draw PROFILE [options] -p <output> -gt <genometable> -i <ChIP>,<Input>,<name> [-i <ChIP>,<Input>,<name> ...]
       <output>: Name of output files
       <genometable>: Tab-delimited file describing the name and length of each chromosome
       -i: specify ChIP data, Input data and name of ChIP sample (separated by ',', <Input> and <name> can be omitted)

Options:
       (以下オプションが続く)

まずはデフォルトで実行してみましょう。コマンドは以下のようになります。

drompa_draw PROFILE \
-i parse2wigdir/H3K4me3,parse2wigdir/Input,H3K4me3 \
-i parse2wigdir/H3K27me3,parse2wigdir/Input,H3K27me3 \
-i parse2wigdir/H3K36me3,parse2wigdir/Input,H3K36me3 \
-p profile.K562 -gene refFlat.txt -gt genometable.txt 

このコマンドにより以下のprofile.K562.pdfが生成されます。

f:id:rnakato:20180524223541p:plain:w350

このプロファイルは遺伝子の転写開始点(TSS)周辺5kbp (±2.5kbp) の平均マップリード数を可視化したものです。ラインが各サンプルの平均リード分布、陰になっている領域が95%信頼区間を表しています。
この図を見ると、H3K4me3はTSS周辺にリードの濃縮がみられる一方、H3K27me3とH3K36me3は濃縮がみられないことがわかります。

オプション

-cwオプションをつけると可視化する周辺長の長さを変えることができます。±5kbpのプロファイルを生成するには以下のようにします。

drompa_draw PROFILE \
-i parse2wigdir/H3K4me3,parse2wigdir/Input,H3K4me3 \
-i parse2wigdir/H3K27me3,parse2wigdir/Input,H3K27me3 \
-i parse2wigdir/H3K36me3,parse2wigdir/Input,H3K36me3 \
-p profile.K562 -gene refFlat.txt -gt genometable.txt -cw 5000

f:id:rnakato:20180524223642p:plain:w350

デフォルトではChIPサンプルのリードの平均値を描画しますが、-stype 1 をつけるとChIP/InputのEnrichmentの平均値を描画します。

drompa_draw PROFILE \
-i parse2wigdir/H3K4me3,parse2wigdir/Input,H3K4me3 \
-i parse2wigdir/H3K27me3,parse2wigdir/Input,H3K27me3 \
-i parse2wigdir/H3K36me3,parse2wigdir/Input,H3K36me3 \
-p profile.K562 -gene refFlat.txt -gt genometable.txt -stype 1

f:id:rnakato:20180524223834p:plain:w350

この例だとChIP readの時との差がわかりにくいですが、broad markや酵母の複製解析などでは差が顕著になることもあります。
(これは各TSSのChIP/Input の平均値であり、ChIP平均/Input平均のEnrichmentではないことに注意。)

遺伝子全長プロファイル

次に、遺伝子全長を100分割したプロファイルを作ってみましょう。 -ptype 3をつけて実行します。

drompa_draw PROFILE -ptype 3 \
-i parse2wigdir/H3K4me3,parse2wigdir/Input,H3K4me3 \
-i parse2wigdir/H3K27me3,parse2wigdir/Input,H3K27me3 \
-i parse2wigdir/H3K36me3,parse2wigdir/Input,H3K36me3 \
-p profile.K562.ptype3 -gene refFlat.txt -gt genometable.txt 

f:id:rnakato:20180519134337p:plain:w400

このプロファイルは各遺伝子について長さを一定に正規化したうえで遺伝子長を100分割し、遺伝子内(gene body)でどのように濃縮しているかを見たものです。 x軸の0がTSS、100がTES(転写終了点)で、その両側に遺伝子長と同じ長さだけ伸長した領域について計測します。 (従って遺伝子ごとに計測している領域の幅が異なることに注意)
100分割した1つのウィンドウがwigファイルの1bin (100bp)より大きい場合、その領域内の全binの平均値を出力します。

TSS周辺の可視化ではほとんど平らだったH3K36me3ですが、遺伝子全体で見るとgene bodyに広くゆるやかな濃縮が広がっていることがわかります。H3K36me3は転写中のgene bodyに入る修飾ですので、これは妥当な結果と言えます。
H3K27me3の方はTSS周辺に微小な濃縮があるようですが、こちらは意味があるかどうか微妙なところです。

領域内でのリード数正規化

3つのサンプル間で縦軸がずれている(平均リード値が異なる)のが気になる場合は、-ntype 1オプションを追加すると、観測している領域内の総リード数でサンプルごとに正規化します。

drompa_draw PROFILE -ptype 3 -ntype 1 \
-i parse2wigdir/H3K4me3,parse2wigdir/Input,H3K4me3 \
-i parse2wigdir/H3K27me3,parse2wigdir/Input,H3K27me3 \
-i parse2wigdir/H3K36me3,parse2wigdir/Input,H3K36me3 \
-p profile.K562.ptype3 -gene refFlat.txt -gt genometable.txt 

f:id:rnakato:20180519151326p:plain:w350

完全ではありませんが、遺伝子外領域の高さがだいたい同じになりました。

まとめ

遺伝子全長の可視化は、H3K36me3のようにgene bodyに入るヒストン修飾や、RNA pol2のようにgene bodyに広く分布する(伸長する)タンパク質のChIP-seqを可視化するのに有効です。 今回はTSS周辺と遺伝子全長について可視化しましたが、-ptype 2でTES周辺、-ptype 4で指定したピークリスト周辺の濃縮を可視化することができます。

また、同時に生成されるRスクリプトprofile.K562.R を編集の後以下のように実行することで、ラベルやy軸のスケールを変更したpdfファイルを再生成することも可能です。

 R --vanilla < profile.K562.R

プロファイルの場合全体の何割で濃縮が入っているかはわかりませんが、信頼区間の幅を見ることである程度推測することが可能です。

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

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

今回は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でそれぞれ出力される図、スタッツのファイル名を指定します。 --ignoreDuplicatesPCR 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の最も大きな問題は、スコアがリード数に依存してしまうという点です。たくさんのサンプルをこの手法で比較する場合には、リード数を揃えてやる必要があります。

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

これは前回の記事の続きですので、未読の方はそちらから読んでください。

Cross-correlation profile

Cross-correlation profile (以下CCP)はENCODEのグループによって提案されたS/N比計測手法です*1。 CCPは、順鎖ー逆鎖間のマップリード分布の「ずれ」を利用します。

f:id:rnakato:20180324131541p:plain:w400

この図で示すように、シーケンサはChIPで得られたDNA断片の片端を固定長で読むため、順鎖と逆鎖にマップされるリード分布にはDNA断片長(約150塩基対)ぶんだけ「ずれ」が発生します。得られるピークも順鎖と逆鎖でずれが起きますので、全てのピーク抽出ツールはこのずれを内部で補正しています。

ということは、この「ずれ」の長さだけ順鎖と逆鎖をずらしてやった時にリード分布が重なるはずです。やってみましょう。順鎖と逆鎖をずらしながら、互いのリード分布の相関を取ったものが下図です。

f:id:rnakato:20180322195938p:plain:w400

これが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を見てみましょう。

f:id:rnakato:20180323121908p:plain:w400

ちゃんとDNA断片長のところが最大になっています。リード長部分のスパイクはほとんど見えないので、これはよいサンプルですね。
下に表示されている"strand-shift (130)"が推定されたDNA断片長です。約130塩基対ということなので、やや短めです。
NSC, RSCS/N比のスコアで、NSCはバックグラウンドに対するDNA断片長(赤点線)部のピーク高さ、RSCはリード長(青点線)部に対するDNA断片長(赤点線)部のピーク高さを表します。 QtagNSCとRSCの値を元に判定されるサンプルクオリティの(S/N比)のレベルを表しています。-2,-1,0,1,2の5つの値を持ち、2が最高、-2が最低です。このサンプルは2なので、最高レベルです。

続いてInputサンプルの方を見てみましょう。

f:id:rnakato:20180323122901p:plain:w400

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を解説しようと思います。長くなりますね…。

*1:Landt, S. G. et al. (2012). ChIP-seq guidelines and practices of the ENCODE and modENCODE consortia. Genome Res, 22(9), 1813–31.

*2:Carroll, T. S. et al. (2014). Impact of artifact removal on ChIP quality metrics in ChIP-seq and ChIP-exo data. Front Genet, 5, 75.

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

私が開発したChIP-seqの品質評価ツール"SSP"の論文がBioinformatics誌にアクセプトされました。

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

本論文はbioRxivでも無料公開されています。

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

SSPは、ChIP-seqサンプルのS/N比の計測の他、得られたピークの信頼性、ピークのモード(sharp/broad)の評価などが可能です。

論文出版記念ということで、このSSPの利用法とメリットについて解説したいと思いますが、そのためにはまず、ChIP-seqサンプルのS/N比計測のための既存手法について解説する必要がありますので、今回はそれについて述べます。

S/N比

S/N比(signal-to-noise ratio)は抗体の力価を表しますが、ChIP-seq解析においては「ゲノム全体でどの程度ピークが得られているか」を表します。平たく言えば、そのサンプルから得られたピークの総数及び強度(ピークの高さ)の程度です。
例えば同一細胞種・同一抗体で生産したreplicate間でもテクニカルな要因により得られるピーク数が大きくばらつく場合がありますが、この場合S/N比が異なるサンプルと表現します。 当然ながらChIPサンプルはたくさんピークが取れていて欲しいのでS/N比は高い方が良く、Inputサンプルはピークが多い=ノイズが多いことを意味しますのでS/N比が低い方が望ましいです。

このS/N比の計測手法については、簡単なようで色々と問題点があり、全てのサンプルに対して適用可能な方法が今まで存在しませんでした。SSPはこれらの課題をクリアすべく開発したものです。以下、詳細です。

S/N比の計測手法

最も簡単なS/N比の計測手法は「ピークの数をカウントする」というものですが、これは同一抗体を使った場合のように、ピークの形がサンプル間で共通であるという前提が必要になります。例えばH3K4me3とH3K36me3のようにsharpなピークとbroadなピークを持つサンプルは、単純なピーク数では比較できません。また、broadなピークを持つサンプルでも「broadさ」はかなり違いますので、broad-broad間でもやはり比較できません。

そこで、FRiP (fraction of reads in peaks) というスコアが提案されました*1。 これもやはり単純で、「得られたピーク領域内にマップされたリード数の全マップリード数に対する割合」を表します。 マップリード数が100万で、ピーク領域内にマップされたリード数が20万であれば、FRiPの値は0.2になります。スコアが高いほどS/N比が高いことになります。
この方法だとサンプル間のピーク形状の違いを吸収することができますが、ピークの数と強度は区別できません。 つまり、「ピーク数は少ないが強度が強いサンプル」と「ピーク数は多いが強度は弱いサンプル」が同じFRiPスコアを持つことはあり得ます。

更に、ピーク数にしろFRiPにしろ、共通の問題点があります。それは、得られるスコアの値がピーク抽出ツールやマップリード数に依存するということです。
これらのスコアは実際に得られたピーク数に基づいて計算しますが、ピーク数はピーク抽出の方法を変えれば当然変化します。例えばInputサンプルを用いた場合と用いない場合では、得られるピークセットは大きく変わるでしょう。異なるInputサンプルを使った場合でもピーク数は変わるかもしれません。「ざっくり知ることができれば十分」という場合は別によいのですが、厳密に評価したい場合には問題です。
また、同一のサンプルでもリード数を増やせば増やすほど統計的有意になる領域が増えるため、得られるピーク数もそれに伴って増加します。従って、FRiPの値をサンプル間で比較する際にはサンプル間であらかじめリード数を統一する(ダウンサンプリングする)必要があります。ダウンサンプリングは時間的にも手間的にも面倒で、サンプル数が多い場合には極めて大変です。が、逆にそういった大規模な解析(データベース作成など)ほどS/N比の計測など品質評価を行いたいわけです。可能であれば、ツールやリード数に依存しないS/N比評価を可能にしたいのです。

そこでCross-correlation profileが登場するのですが、長くなってしまいましたので、日を改めて書きたいと思います。 ちなみに今日書いたような品質評価の話は私が以前執筆したレビュー論文*2によくまとまっていますので(自画自賛)より詳細が知りたいというかたは参照されてください。

*1:Landt, S. G. et al. (2012). ChIP-seq guidelines and practices of the ENCODE and modENCODE consortia. Genome Res, 22(9), 1813–31.

*2:Nakato, R. and Shirahige, K. (2017). Recent advances in ChIP-seq analysis: from quality management to whole-genome annotation. Brief Bioinform, 18(2), 279–290.

DROMPA3: その8 GVコマンドでのマクロな可視化

今回は、全染色体を1行でマクロに可視化するGVコマンドを使います。なおGVはGlobal viewの略です。

parse2wig

今回はROADMAP web portalからダウンロードしたK562細胞のヒストン修飾データ一式を使います。 以下のコマンドでtagAlignファイルをダウンロードします。 (ダウンロードに時間がかかる場合は適宜サンプル数を減らして大丈夫です。)

$ wget http://egg2.wustl.edu/roadmap/data/byFileType/alignments/consolidated/E123-H3K4me3.tagAlign.gz
$ wget http://egg2.wustl.edu/roadmap/data/byFileType/alignments/consolidated/E123-H3K9me3.tagAlign.gz
$ wget http://egg2.wustl.edu/roadmap/data/byFileType/alignments/consolidated/E123-H3K27ac.tagAlign.gz
$ wget http://egg2.wustl.edu/roadmap/data/byFileType/alignments/consolidated/E123-H3K27me3.tagAlign.gz
$ wget http://egg2.wustl.edu/roadmap/data/byFileType/alignments/consolidated/E123-H3K36me3.tagAlign.gz
$ wget http://egg2.wustl.edu/roadmap/data/byFileType/alignments/consolidated/E123-Input.tagAlign.gz

ダウンロードしたファイルからparse2wigでWigデータを生成します。 ROADMAPのデータはgenome build hg19です。

for name in H3K27ac H3K27me3 H3K36me3 H3K4me3 H3K9me3 Input
do
    parse2wig -i E123-$name.tagAlign.gz -o E123-$name -gt genometable.txt -f TAGALIGN -binsize 100000 
done

# 上のfor文は以下のコマンドを実行したのと同じ
$ parse2wig -i E123-H3K27ac.tagAlign.gz -o E123-H3K27ac -gt genometable.txt -f TAGALIGN -binsize 100000
$ parse2wig -i E123-H3K27me3.tagAlign.gz -o E123-H3K27me3 -gt genometable.txt -f TAGALIGN -binsize 100000
$ parse2wig -i E123-H3K36me3.tagAlign.gz -o E123-H3K36me3 -gt genometable.txt -f TAGALIGN -binsize 100000
$ parse2wig -i E123-H3K4me3.tagAlign.gz -o E123-H3K4me3 -gt genometable.txt -f TAGALIGN -binsize 100000
$ parse2wig -i E123-H3K9me3.tagAlign.gz -o E123-H3K9me3 -gt genometable.txt -f TAGALIGN -binsize 100000
$ parse2wig -i E123-Input.tagAlign.gz -o E123-Input -gt genometable.txt -f TAGALIGN -binsize 100000

入力がtagAlign形式データのため、オプションに -f TAGALIGNを指定しています。 また、GVコマンドはデフォルトbinsizeが100kbpのため、-binsize 100000 を指定し、100kbpのWigファイルを生成します。

DROMPAの実行

それでは可視化してみましょう。GVコマンドもPC_SHARPコマンドと指定方法はほぼ同様です。

$ drompa_draw GV \
$ -i parse2wigdir/E123-H3K4me3,parse2wigdir/E123-Input,H3K4me3 \
$ -i parse2wigdir/E123-H3K27ac,parse2wigdir/E123-Input,H3K27ac \
$ -i parse2wigdir/E123-H3K27me3,parse2wigdir/E123-Input,H3K27me3 \
$ -i parse2wigdir/E123-H3K36me3,parse2wigdir/E123-Input,H3K36me3 \
$ -i parse2wigdir/E123-H3K9me3,parse2wigdir/E123-Input,H3K9me3 \
$ -p K562 -gt genometable.txt

以下のような図が生成されます。 f:id:rnakato:20180216144039p:plain

これはChIP/Input比を100kbpの解像度で可視化したもので、ChIP/Input >1.0の領域はオレンジ、そうでない領域は灰色で表示されます。
ヘテロクロマチンマーカであるH3K9me3のS/N比は非常に弱いため、通常のピーク抽出ではほとんどピークが抽出されないこともしばしばありますが、このようなマクロな可視化で見ると、他のアクティブマーカであるヒストン修飾と排他的に局在するさまが視覚的によくわかります。染色体の免疫染色のようなイメージですね。

オプションをつければPC_SHARPと同様にリード数そのものを可視化することももちろん可能です。(以下の記事参照)

DROMPA3: その6 ChIP/Input ratio 及び p値の可視化 - Palmsonntagmorgen

GC含量と遺伝子密度のグラフを追加

上の図だと少し殺風景なので、GC含量と遺伝子密度のグラフを追加してみましょう。 hg19のGC含量と遺伝子密度のデータはDROMPAのweb siteからダウンロード可能です。

$ wget http://www.iam.u-tokyo.ac.jp/chromosomeinformatics/rnakato/drompa/annotations/GC_hg19.tar.gz
$ wget http://www.iam.u-tokyo.ac.jp/chromosomeinformatics/rnakato/drompa/annotations/gene_density_hg19.zip
$ tar zxvf GC_hg19.tar.gz
$ unzip gene_density_hg19.zip

GC_hg19は500kbpごとのゲノム配列のGC含量、gene_density_hg19はその領域に含まれる遺伝子数です。 これらをプロットするには以下のように指定します。

$ drompa_draw GV \
$ -i parse2wigdir/E123-H3K4me3,parse2wigdir/E123-Input,H3K4me3 \
$ -i parse2wigdir/E123-H3K27ac,parse2wigdir/E123-Input,H3K27ac \
$ -i parse2wigdir/E123-H3K27me3,parse2wigdir/E123-Input,H3K27me3 \
$ -i parse2wigdir/E123-H3K36me3,parse2wigdir/E123-Input,H3K36me3 \
$ -i parse2wigdir/E123-H3K9me3,parse2wigdir/E123-Input,H3K9me3 \
$ -p K562_withGCandGene -gt /home/Database/UCSC/hg19/genome_table \
$ -GC GC_hg19 -gcsize 500000 \
$ -GD genedensity -gdsize 500000

-GC, -GD オプションでそれぞれのデータ(ディレクトリ)を指定します。-gcsize, -gdsizeで、ウィンドウサイズを指定します。ここでは500kbpです。

結果は以下のようになります。

f:id:rnakato:20180216144049p:plain

H3K9me3が濃縮している領域は遺伝子数の少ないgene desertのような領域が多いことがざっと観察できます。 GC含量は、例えばリードの濃縮が極端にGC richに偏っているかどうか調べたい場合などに便利です。

ちなみに、EnsemblUCSC genome browserで使われている↓のような簡単な染色体図も同時に可視化できるようにしたいのですが、著作権フリーな染色体模式図の画像ってどこかにあるんでしょうか・・・ご存知の方がおられましたら教えてください。

f:id:rnakato:20180216145832p:plain

補足:GC含量データ、遺伝子密度データの作成

GC含量、遺伝子密度のデータは自分で作成することもできます。ここではヒトゲノムhg19を例に解説しますが、他のゲノムの場合には適宜読み替えてください。

GC含量はDROMPA3のscriptフォルダに含まれる GCcount.plを使います。ここではGC_hg19フォルダを新たに作り、その中にデータを生成していくことにします。

$ mkdir GC_hg19  # フォルダ作成
$ winsize=500000  # ウィンドウサイズを500kbpに指定
$ for i in $(seq 1 22) X Y M # 1~22番染色体とX, Y染色体、ミトコンドリアを指定。
$ do
$    GCcount.pl chr$i.fa $winsize > GC_hg19/chr${i}_bs$winsize
$ done

GCcount.pl は単一のfastaファイルを入力とするので、染色体ごとに適用します。 出力ファイル名は<染色体名>_bs<ウィンドウサイズ>である必要があり、染色体名はgenometable.txtと一致している必要があります。
このように作成し、-GC GC_hg19と指定すれば可視化可能です。任意のウィンドウサイズで可視化可能で、PC_SHARPコマンドでも可視化できます。

遺伝子密度データはDROMPA3のscriptフォルダに含まれる makegenedensity.pl を使います。 実行にはgenome tableと、refFlat形式の遺伝子ファイルが必要です。持っていない場合はUCSCからダウンロードしましょう。

# 遺伝子ファイルのダウンロード 
$ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/database/refFlat.txt.gz
$ gunzip refFlat.txt.gz  # 解凍
# 遺伝子密度データ作成
$ winsize=500000  # ウィンドウサイズを500kbpに指定
$ makegenedensity.pl genometable.txt refFlat.txt $winsize

入力にするrefFlatファイルをあらかじめ編集しておくことで、例えばコーディング遺伝子のみ・ノンコーディング遺伝子のみカウントするような使い方も可能です。