Palmsonntagmorgen

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

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

NGS界隈におけるプログラミング言語の競争について – 極めて主観的な見地から(1/30・31追記)

タイトルはこの有名な記事からもらいました。
自分の学生時代に講義で学んだプログラミング言語はCとPerl(とjavaPostgreSQL)でしたが、状況はずいぶんと変わってきました。じゃあどういう時代になったの?というのを、自分が今いるNGS解析のフィールドから見える景色として記録しておきます。 自分自身がこの分野の最先端にいるという気はしないので見当違いもあるかと思いますが、酒の肴と思ってご笑覧ください。 なお、これは論文として発表されるオープンソースツールについての記述が主で、商用ツールに関しては触る機会がほとんどないのでよくわかりません。

(1/31追記)こういうエントリを書くと世の中の他の言語が気になりだしますね。。これからも増える可能性があるので、その他の新言語については下記に順次追加するようにします。

C

新しく論文で発表された解析ツールがゴリゴリのCで書かれているなんてことは流石になくなってきた。 とはいえ、高速化・省メモリ化の点では最強なので、ソースは巨大にならないが入力データが巨大であるような分野では選ばれるかもしれない。流れとしてはC++に移行しているんだけど、C++よりCのが速いよね??昔から存在するCで書かれた解析ツール(bedtoolsとか)は最近CとC++のハイブリッド的な感じで更新しているようだ。
個人的にはCを書いている時が一番楽しい。既存ツールをコンパイルする時にgccが動いていると心があたたまる。

C++

Bowtie, Hisat, Salmon, STARなど、高速性が何よりも重要であるゲノムマッピングの分野ではメインになっている。Cよりも柔軟にソースが書けるのでコードがすっきりする。ライブラリもCより充実してて嬉しい。コーディング的な意味で独創的なアイデアを詰め込みやすい。俺が学生の時のC++は単なるオブジェクト指向Cだったんだけど、テンプレートプログラミングが登場してからは言語そのものがクラスチェンジした感あり。抽象的な構造になるため、エレガントなプログラムソースはもはや知恵の輪のような、綺麗だけど読めねー的な要素を醸しだす。C++で書いた解析ツールを発表するのはプログラミング能力に自信がある人達という印象。コンパイルでもCMakeを使いこなす。 C++17, C++20とどんどん進化していっているんだけど、そういう最先端のC++で書くようなツールは今後出てくるのだろうか?ちなみにDROMPAはCですがそのうちC++にクラスチェンジします。

R

統計処理で最も威力を発揮するので、複雑な確率モデルを必要とするRNA-seqの発現解析でよく利用されている。モデリングや検定をコマンド一行でできるのずるい。中身理解してなくても複雑なモデリングできちゃうじゃん。Bioconductorが非常に発展していて、ガラパゴス化しがちだったツール同士の統合がすごく楽になった。今のBioinformaticsの発展を支える一助になっていると思う。一方やたら低速高メモリ消費なのででかい行列とかを食わせるとハングする。(それなのに全ゲノムのHI-Cマトリクスを食わせるツールとか正気じゃないぜ。)マルチスレッディングも苦手。という訳で遅い。大規模解析にはあまり向いてない。ChIP-seqのようにノイズがまじりやすいデータにもあまり向いてないと思う(場合によるけど)。CやPythonで書かれたツールを全部ラップしてRツールにしてしまう癖がある。Rでディープラーニングができます!って銘打ったツールを使ってみたら内部でKerasを起動してるだけでワロタ。

Perl

Perlは死んだとか言われてるけど、死んでねーし。 流石に今からプログラミング学びたいって人に薦める気はないけど、既にPerlで書かれてうまく動いているプログラムをPythonに書き直す必要はないでしょ?awkより高機能なパーサだと思えば便利。ワンライナーも書けるし、文字列のハッシュも便利。NGS解析ツールの分野でもPerlで書かれた現役ツールはけっこうあったりする。HOMERが最たる例。.plの拡張子消してバイナリファイルっぽくしてるけどこっそりPerlってツールはけっこうある、RSEMとか。マイナーどころではSISSRとか。でも正直、新規ツールの論文見てソース参照したらPerlスクリプトだったりすると、残念な印象を受けることは否めない。 そういえばBioPerlってのが昔あったんだけど、あれを実際に自分の研究に利用してた/してる人ってどのくらいいたんだろう。。。

Python

Pythonの歴史は古いが爆発的に流行りだしたのはPandasが登場したからだと言われている。Pandasはテキストパーサとして非常に便利。流行りのGPUコンピューティングディープラーニングするならPython一択みたいなところがあるので、その界隈の人たちは皆Pythonユーザである。Cythonだから中身はCだけどな。ツールのインストールが初心者にとって最初のつまずきポイントであることは異論が無いだろうが、pip, condaでインストールが楽になったことも流行る一因だろう。あとは可読性が高いから、とかあるんだけど、実際のところどうかな。 プログラミングの知識自体はそこまで無くともコピペ改変でプログラムが書けてしまうので、素人っぽいソースがたくさんある。 なおPython2系は2020年でサポートが切れるのだが、NGS解析ツールはPython2ばかり。大丈夫なのだろうか。 ごくたまにPython3で作られたプログラムも見かける。そういう人たちは実行にSnakemakeを使うんだけど、なんかわかりにくいのでやめて欲しい。なお作図能力ではRに一歩劣る。

Ruby

使ったことなし。Perlっぽい感じ?Railsとして生き残っている?

Java

GUIを持つブラウザなどを除くとあまりJavaツールを見かける機会はないのだが、一方でChromHMM, Juicer, GATKなど、各分野でのキラーアプリに近いツールがJavaを使って実装されているのは注目に値する。理由はよくわからないが、内製のライブラリとかが研究室(大学?)にあるのかな?本来はインタラクティブな可視化ブラウザなどを実装する時に威力を発揮する言語。ちなみにJavaScriptとは別物です。

MATLAB

機械学習アルゴリズムを実装するのにとても便利なのでよく使われるが、有償なので、いずれはPythonとRに取って代わられるだろう。と言われているのだが、おそらくNGSの分野ではずっと現役だと思う。何故ならMATLABで書かれた新規ツールは今でもよく見かけるし、MATLABがメイン言語の研究室はおそらく相当数ある。これまでの資産がMATLABで、新人の教育フローもMATLABでうまく回っている研究室がそう簡単にMATLABを手放すとは思えない。Python2もだけど、研究室のメイン言語を切り替えるのってたぶん世の中の流れより5年は遅い。なので、既存ツールを全部使えるようにしておきたいのならこれからもMATLABは捨てられないことになるけど、うーん、今から新人に学ばせるのもなあ。

Octave

かのNg大先生おすすめの言語だけど、あんまり見ることはない。使ったこともないです。MATLABと同じなのかな?

(1/30追記)Rust

Rust使ってる人っていますか?僕はこの記事書いた時点で、この言語を知りませんでした(名前は見かけていたけど、プログラミング言語だと認識していなかった)。ポインタの苦労から解放されたC++という話を聞きましたが、どうなんでしょう。使ってみましょうかね。DROMPAくらい大きなプログラムには難しいでしょうが、小さいツールならばできそうです。しかし次々新しい言語ってできるもんですね…

その他の言語 (1/31追記)

Julia、lua, Go, Haskell, Nimなど。注目の言語だけど今のところこれらを使って実装されている解析ツールは見たことがない。でもin-houseで使っている人はおそらくいるのだろう。基本的な統計処理、検定などの実装には使えそうだけど、例えばNGSのBAMファイルを処理したりというようなことは難しい(あるいは低速)と思うので、他の言語と併用で使うか、予備実験用に使うなどの用途になるだろう。Juliaはユーザとライブラリが増えればPythonのように伸びる可能性はある。わからんけど。

結論

NGS分野でのプログラミング言語の流行は世の中の流れと微妙に違っていて、かつ混沌としている。例えばBioconductorのプログラムを全部Anacondaに移植しようとしてる猛者なんかもいるんだけど、結局エラーとの闘いになるし、Anacondaも先行き不透明だったりするので、自分の必要な分野の動向を見つつ、適当にいくつか使えるようにしておくというのが良い方法なのだろう。論文出したら後は知らん、というスタンスならどれで書いてもいいんだけど、実際に使ってもらうために何年もメンテナンス・アップデートするつもりであれば、どの言語で書くか(あるいは移植するか)は重要な問題になるだろう。

BEDtoolsワンライナー覚書

BEDtoolsの作者は開発熱心なので、できることがどんどん増えているような気がします。
手元のバージョンはv2.27.1です。

前準備

多くのコマンドはsorted BEDを要求しますので、事前に以下のコマンドで全てのBEDをソートしておくとストレスがないかと思います。

 $ sort -k1,1 -k2,2n in.bed > in.sorted.bed

更に、いくつかのコマンドにはGenome tableファイルが必要になります。 Genome tableの作成は以下の記事を参照してください。

genome tableを作成する - Palmsonntagmorgen

a.bedに含まれる領域のうち、b.bedと重なる領域"のみ"を出力

 $ bedtools intersect -a a.bed -b b.bed > a_and_b.bed

-bには複数のBEDファイルを指定することが可能です。

 $ bedtools intersect -a a.bed  -b b1.bed b2.bed b3.bed > a_and_b.bed

a.bedに含まれる領域のうち、b.bedと重なる領域を出力

$ bedtools intersect -wa -a a.bed -b b.bed > a_and_b.bed

intersectで出力される領域については、以下に図解があります。
intersect — bedtools 2.27.0 documentation

a.bedに含まれる領域のうち、b.bedと重ならない領域を出力

$ bedtools intersect -v -a a.bed -b b.bed > a_not_b.bed

a.bedの各領域がb.bedにどの程度カバーされるかを調べる

$ bedtools coverage -a a.bed -b b.bed > a.coverage.bed

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

chr1  0   100  3  30  100 0.3000000
chr1  100 200  1  100 100 1.0000000
chr2  0   100  0  0   100 0.0000000
...

1~3列目はa.bedと同一です。4列目はそのサイトと重なるb.bed内のサイトの数、5列目はサイト長に対するb.bedがoverlapしている領域の長さ、6列目はサイト長、7列目は割合 (5列目/6列目)です。

BEDに含まれる領域を両側に100bp伸ばす

 $ bedtools slop -i in.bed -g genome_table -b 100 > in.extend.bed

片側のみに伸長したい場合は-l, -rオプションを使います。

BEDに含まれる領域のうちoverlapするものをマージ

 $ bedtools merge -i in.bed > in.merged.bed

-dオプションでoverlapの基準をゆるめることができます。

BEDに含まれる領域を両側に1000bp伸ばし、重なったサイトはマージ

$ bedtools slop -i in.bed -g genome_table -b 1000 | bedtools merge -i - > in.extend.bed

3つのBEDファイルを結合し、重なるサイトはマージ

catするとsortされていないBEDが出力になるので、間にsortをかませています。

$ cat a.bed b.bed c.bed | sort -k1,1 -k2,2n | bedtools merge -i - > merged.bed

BEDに含まれて’いない’ゲノム領域を出力

$ bedtools complement -i in.bed -g genome_table > notin.bed

全ゲノム領域について、 BEDに含まれた領域にどの程度含まれているかをヒストグラムで出力

3つ以上のピークリストに対して、どの程度ピーク同士が重なっているかを調べる時などに使います。 BAMを入力にしてリードをカウントすることもできますが、遅いです。

$ bedtools genomecov -i in.bed -g genome_table > in.genomecov.bed

出力については以下の図解がわかりやすいです。

genomecov — bedtools 2.27.0 documentation

ゲノム内で100bpの領域を10000個ランダムに生成

ゲノムにはギャップ領域などがありますので、実際にはそれらを除外した領域を与えることが望ましいでしょう。

$ bedtools random -l 100 -n 10000 -g genome_table > random.bed

ランダムのため得られる結果は毎回変わりますが、-seedオプションを指定することで再現可能になります。

in.bedに含まれる領域をゲノム内でシャッフルして出力

randomではサイト長が固定ですが、こちらは実際に得られたピークデータを入力にするぶん、より実際に近いと思います。
Permutation testなどで使えそうですが、randomと同様、与えるゲノム領域を考慮する必要があります。

$ bedtools shuffle -i in.bed -g genome_table

シャッフルを同一染色体内に限定

$ bedtools shuffle -i in.bed -g genome_table -chrom

シャッフル後のサイトをinclude.bedに含まれる領域のみに限定

$ bedtools shuffle -i in.bed -g genome_table -incl include.bed

シャッフル後のサイトがexclude.bedに含まれる領域に含まれないようにする

$ bedtools shuffle -i in.bed -g genome_table -excl exclude.bed

シャッフル後のサイトが互いに重ならないようにする

$ bedtools shuffle -i in.bed -g genome_table -noOverlapping

BEDで指定されている領域の塩基配列を取得

モチーフ解析の時に必要になります。

$ bedtools getfasta -fi genome.fa -bed in.bed > in.fa

a.bedとb.bedの重なりをJaccard indexで計算

Jaccard indexとは、平たく言えばベン図の重なり部分の大きさを0~1の間のスコアで評価するものです。 便利そうで、あまり使う機会がありません。

$ bedtools jaccard -a a.bed -b b.bed

詳細は以下を参照してください。

jaccard — bedtools 2.27.0 documentation

複数のBAMファイルに対して、target_regions.bed に含まれる領域のマップリード数をカウント

 $ bedtools multicov -bams aln1.bam aln2.bam aln3.bam -bed in.bed > in.count.bed

parallelコマンドがインストールされていれば、以下のコマンドでも可能です。

 $ find *bam | parallel 'bedtools coverage -hist -abam {} -b target_regions.bed | grep ^all > {}.hist.all.txt'

これは以下のコマンドを実行したのと同じ意味になります。

 $ bedtools coverage -hist -abam samp.01.bam -b target_regions.bed | grep ^all > samp.01.bam.hist.all.txt
 $ bedtools coverage -hist -abam samp.02.bam -b target_regions.bed | grep ^all > samp.02.bam.hist.all.txt
 $ bedtools coverage -hist -abam samp.03.bam -b target_regions.bed | grep ^all > samp.03.bam.hist.all.txt
 $ bedtools coverage -hist -abam samp.05.bam -b target_regions.bed | grep ^all > samp.05.bam.hist.all.txt
 $ bedtools coverage -hist -abam samp.06.bam -b target_regions.bed | grep ^all > samp.06.bam.hist.all.txt
 $ bedtools coverage -hist -abam samp.10.bam -b target_regions.bed | grep ^all > samp.10.bam.hist.all.txt

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

DROMPA3: その4 マップリード分布の可視化その1 では以下のコマンドを実行しました。

$ drompa_draw PC_SHARP \
$ -i parse2wigdir/H3K4me3,parse2wigdir/Input,H3K4me3,,,200 \
$ -i parse2wigdir/H3K27me3,parse2wigdir/Input,H3K27me3,,,10 \
$ -i parse2wigdir/H3K36me3,parse2wigdir/Input,H3K36me3,,,10 \
$ -p K562 -gene refFlat.txt -gt genometable.txt -chr 1 -lpp 2

今回は-iで指定するデータの詳細について解説します。

-iの内容

-i で指定するのは、ChIPサンプルと対応するInputサンプルのペアのファイル名、及びそのペアに対するオプションをカンマ区切りで指定します。 具体的には以下の並びでDROMPAに考慮されます。1のChIPサンプル名以外は省略可ですが、カンマは省略できません。

  1. ChIPサンプルファイル名(必須)
  2. Inputサンプルファイル名
  3. ChIPサンプルラベル(図中で表記される名前)
  4. Peak BEDファイル(図中でハイライトされる領域)
  5. ビンサイズ
  6. y軸スケール(read行)
  7. y軸スケール(ChIP/Input ratio行)
  8. y軸スケール(p値行)

つまり以下の指定は、

$ -i parse2wigdir/H3K4me3,parse2wigdir/Input,H3K4me3,,,200 

ChIP・Inputサンプルファイル名はそれぞれparse2wigdir/H3K4me3parse2wigdir/Inputであり、ラベルはH3K4me3、read行のy軸スケールは200で、他の値についてはデフォルト値または全体のオプションで指定した値を用いるということになります。

それぞれの詳細

Peakファイルを指定した場合、BEDファイルで指定された領域を赤色でハイライトするようになります。指定しない場合、DROMPAが内部でピークコールを行います。これにより、例えば複数のピーク抽出ツールで得られたピーク領域の比較などができるようになります。

ビンサイズの指定は、例えばsharp peakは100bpだが、broad histone markについては1kbpで可視化したいというような場合に有用です。対応するInputファイルもそのbinsizeのwigデータが必要であることに注意してください。

フルで指定すると、以下のような感じになります。

$ -i parse2wigdir/H3K4me3,parse2wigdir/Input,H3K4me3,H3K4me3_peak.bed,1000,200,2,10