Palmsonntagmorgen

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

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 

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

リード分布の可視化の続きです。 このエントリは↓の記事の続きになりますので、まだ読んでいない方は先にこちらを参照してください。

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

Input readの可視化

前回はChIPサンプルのみを可視化しましたが、Inputサンプルのリード分布を並列して可視化することもできます。 以下のように、-show_itag 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_1 -gene refFlat.txt -gt $gt -chr 1 -lpp 2 -show_itag 1 -rmchr

f:id:rnakato:20171221163620p:plain

各ChIPサンプルの下に、対応するInputサンプルのリード分布が表示されるようになりました。 Inputサンプルのy軸のスケールはChIPサンプルで指定したものに準拠します。

今回の例の場合はInputサンプルは3サンプル共通ですので、別々に表示する必要はありません。 -show_itag 2として、最下段にのみInputサンプルを表示するようにしましょう。

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_2 -gene refFlat.txt -gt $gt -chr 1 -lpp 2 -show_itag 2 -rmchr

f:id:rnakato:20171221163809p:plain

すっきりと表示することができました。なお、Inputサンプルがサンプル毎に異なる場合に-show_itag 2を付加すると、最初に指定されたInputサンプルを代表として可視化します。

ChIP/Input ratio, p値の可視化

DROMPAはリード分布を可視化するだけではなく、ChIP/Input ratio, ピークコールに使った二種類の検定のp値(こちらの記事参照)も可視化することができます。

  • -showratio 1で ChIP/Input ratio,
  • -showpinter 1 でChIP internal p-value,
  • -showpenrich 1 でChIP/Input enrichment p-valueをそれぞれ可視化します。

それではすべてを可視化してみましょう。段数が多くなるので、-lppの値を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_3 -gene refFlat.txt -gt $gt -chr 1 -lpp 1 -show_itag 2 -rmchr \
-showratio 1 -showpinter 1 -showpenrich 1 

f:id:rnakato:20171221164459p:plain

それぞれのサンプルについて、上から ChIP internal、ChIP/Input enrichment p-value、ChIP/Input ratio、ChIP readの順に表示されています。 2種類のp値は-log10(p)の値に変換されています。すなわち、例えばp=0.01の場合は値が2になります。
それぞれのパラメータについて、設定された閾値を上回る領域は赤色、それ以外の領域は灰色で表示されます。上記コマンドではデフォルト設定のため、

  • ChIP internal p < 1e-4
  • ChIP/Input enrichment p < 1e-4
  • ChIP/Input ratio > 0

の領域がそれぞれ赤で表示されます。なお、全ての閾値をクリアした領域、すなわちピーク領域はChIP readのラインが赤色になります*1

このような表示は、適切な閾値の推定に役立ちます。例えばH3K27me3やH3K36me3のサンプルは、濃縮領域と思われる領域が緑色になっており、ピークとして検出されていません。p値のラインを見ると、 ChIP internalの閾値はクリアしているが、ChIP/Input enrichmentのp値が閾値を下回っていることがわかります。従ってこの領域をピークとして検出できるように閾値を緩めたい場合は、ChIP/Input enrichmentの閾値-pthre_enrich)を緩めればよいことになります。 DROMPA PC_SHARPのデフォルトパラメータはH3K4me3のようなstrong sharp peakに最適化されていますので、broad peakやS/N比の弱い抗体などを使った場合は、閾値を緩めた方が良い結果が得られるでしょう。

*1:厳密には、可視化ではq値を考慮していないため、実際に出力されるピークリストと赤色でハイライトされる領域が微妙に食い違う場合があります。また、可視化領域に対してピーク領域が微小すぎる場合、赤色がつぶれて見えなくなる場合があります。

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

今日はシェル変数について。

前回の記事の「複数サンプルの可視化」の項で、以下のコマンドを実行しました。

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

ここで、”parse2wigdir”を6回も記入するのは面倒ですので、これを変数に格納してみましょう。

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

1行目でdirという変数に"parse2wigdir"という文字列を格納しています。 以後、$dirと指定することで"parse2wigdir"を記入したのと同じ意味になります。
(変数を定義する時には$記号がなく、参照する時には$記号が必要であることに注意)
仮にディレクトリ名が変更された場合、最初の例では6箇所修正する必要が出てきますが、このように変数に格納すると修正は1箇所で済みます。 また、「6箇所に同じ名前を利用する」という意味合いがはっきりするという意味でもこちらの方が望ましいです。

もう少し修正してみましょう。以下のようにするとどうでしょうか。

$ dir=parse2wigdir
$ s1="-i $dir/H3K4me3,$dir/Input,H3K4me3"
$ s2="-i $dir/H3K27me3,$dir/Input,H3K27me3"
$ s3="-i $dir/H3K36me3,$dir/Input,H3K36me3"
$ drompa_draw PC_SHARP $s1 $s2 $s3 -p K562 -gene refFlat.txt -gt genometable.txt -chr 1 -lpp 2

3つのサンプルペアをそれぞれs1, s2, s3に格納し、コマンドを実行する時にそれを指定しています。 こうすると以下のようにサンプルを柔軟に選択できるようになります。

$ dir=parse2wigdir
$ s1="-i $dir/H3K4me3,$dir/Input,H3K4me3"
$ s2="-i $dir/H3K27me3,$dir/Input,H3K27me3"
$ s3="-i $dir/H3K36me3,$dir/Input,H3K36me3"
$ drompa_draw PC_SHARP $s1 -p K562_1 -gene refFlat.txt -gt genometable.txt -chr 1 -lpp 2
$ drompa_draw PC_SHARP $s2 -p K562_2 -gene refFlat.txt -gt genometable.txt -chr 1 -lpp 2
$ drompa_draw PC_SHARP $s3 -p K562_3 -gene refFlat.txt -gt genometable.txt -chr 1 -lpp 2
$ drompa_draw PC_SHARP $s1 $s2 $s3 -p K562_all -gene refFlat.txt -gt genometable.txt -chr 1 -lpp 2

3サンプル程度だとここまでする必要はあまりないのですが、数十サンプル同時に扱うような場合には威力を発揮します。

パラメータを毎回指定するのも面倒なので、統一的に使いたいパラメータも変数にしてしまいましょう。

$ dir=parse2wigdir
$ s1="-i $dir/H3K4me3,$dir/Input,H3K4me3"
$ s2="-i $dir/H3K27me3,$dir/Input,H3K27me3"
$ s3="-i $dir/H3K36me3,$dir/Input,H3K36me3"
$ param="-gene refFlat.txt -gt genometable.txt -chr 1 -lpp 2"
$ drompa_draw PC_SHARP $s1 -p K562_1 $param
$ drompa_draw PC_SHARP $s2 -p K562_2 $param
$ drompa_draw PC_SHARP $s3 -p K562_3 $param
$ drompa_draw PC_SHARP $s1 $s2 $s3 -p K562_all $param

ずいぶんすっきりしました。 変数の数を増やしすぎると逆に読みにくくなることもあるので、どこまで変数に格納するかは個人の好みによります。

最後に、20サンプルを同時に可視化する場合を考えてみましょう。 ChIPサンプルはChIP1, ChIP2, ..., ChIP20とナンバリングされているとします。Inputサンプルは共通です。 上の例にならうと、以下のようになります。

$ dir=parse2wigdir
$ s1="-i $dir/ChIP1,$dir/Input,ChIP1"
$ s2="-i $dir/ChIP2,$dir/Input,ChIP2"
...
$ s20="-i $dir/ChIP20,$dir/Input,ChIP20"
$ param="-gene refFlat.txt -gt genometable.txt -chr 1 -lpp 2"
$ drompa_draw PC_SHARP $s1 $s2 ... $s20 -p K562_all $param

サンプルの定義に20行も使ってしまいました。これはあまりスマートではありません。 こんな時はfor文を使うと簡単に書くことができます。何が起きているかわかるでしょうか。

$ dir=parse2wigdir
$ s=""    # 空の変数を定義
$ for i in $(seq 1 20)    # i を1から20まで順にインクリメント
$ do
$     s="$s -i $dir/ChIP$i,$dir/Input,ChIP$i"  # 変数sにサンプルを再帰的に追加
$ done
$ param="-gene refFlat.txt -gt genometable.txt -chr 1 -lpp 2"
$ drompa_draw PC_SHARP $s -p K562_all $param

このようなコマンドのかたまりをテキストファイルとして保存したものをシェルスクリプトと言います。 DROMPAの一番の弱点はコマンドが長くなってしまうことなのですが、このようにするとすっきりと記述することができますので、試してみてください。

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

今回はDROMPA3のメイン機能であるマップリード分布の可視化を紹介します。

インストール

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

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

今回は複数のChIPサンプルを可視化したいので、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

drompa_draw

可視化にはdrompa_drawを使います。 "drompa_draw"とタイプすると、利用可能なコマンド一覧が出力されます。

$ drompa_draw
Usage: drompa_draw [--version] <command>
Command: PC_SHARP    peak-calling (for sharp mode)
         PC_BROAD    peak-calling (for broad mode)
         PC_ENRICH   peak-calling (enrichment ratio)
         GV          global-view visualization
         PD          peak density
       ... 

今回はPC_SHARPを使います。 "drompa_draw PC_SHARP"とコマンド含めてタイプすると、そのコマンドで利用可能なオプション一覧が表示されます。

$ drompa_draw PC_SHARP
Usage: drompa_draw PC_SHARP [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:
       (以下オプションが続く)

1サンプルの可視化

簡単のため、最初は1サンプルだけ使って可視化してみます。

$ drompa_draw PC_SHARP -i parse2wigdir/H3K4me3,parse2wigdir/Input,H3K4me3 \
-p K562_H3K4me3 -gene refFlat.txt -gt genometable.txt 

-i でChIP、Inputサンプルをカンマ区切りで指定するのはdrompa_peakcallと同じです。更にその後にコンマを入れて"H3K4me3"とタイプしていますが、これは出力図に表示されるサンプル名です。 -pで出力ファイル名、-geneで遺伝子アノテーションを指定します。

無事成功すると、以下のように全ゲノムを可視化したpdf(K562_H3K4me3.pdf)と、染色体別のpdfファイル(K562_H3K4me3_chr*.pdf)が作成されるでしょう。

$ ls
K562_H3K4me3.pdf        K562_H3K4me3_chr14.pdf  K562_H3K4me3_chr2.pdf   K562_H3K4me3_chr5.pdf
K562_H3K4me3_chr1.pdf   K562_H3K4me3_chr15.pdf  K562_H3K4me3_chr20.pdf  K562_H3K4me3_chr6.pdf
K562_H3K4me3_chr10.pdf  K562_H3K4me3_chr16.pdf  K562_H3K4me3_chr21.pdf  K562_H3K4me3_chr7.pdf
K562_H3K4me3_chr11.pdf  K562_H3K4me3_chr17.pdf  K562_H3K4me3_chr22.pdf  K562_H3K4me3_chr8.pdf
K562_H3K4me3_chr12.pdf  K562_H3K4me3_chr18.pdf  K562_H3K4me3_chr3.pdf   K562_H3K4me3_chr9.pdf
K562_H3K4me3_chr13.pdf  K562_H3K4me3_chr19.pdf  K562_H3K4me3_chr4.pdf   K562_H3K4me3_chrX.pdf

染色体別ファイルが必要ない時は-rmchrオプションを追加すると全ゲノムpdfのみ生成します。 また、例えば染色体1番のみで良い場合は、以下のように-chr 1と追加するとよいでしょう。

$ drompa_draw PC_SHARP -i parse2wigdir/H3K4me3,parse2wigdir/Input,H3K4me3 \
-p K562_H3K4me3 -gene refFlat.txt -gt genometable.txt -chr 1
...
$ ls
K562_H3K4me3.pdf   K562_H3K4me3_chr1.pdf  

K562_H3K4me3_chr1.pdf を開いてみましょう。

f:id:rnakato:20171212183101p:plain

上部は遺伝子領域です。中央より上にあるものは順鎖、下にあるものは逆鎖です。refFlatデータはisoformを区別しますので、表示される遺伝子もtranscript単位になっています(ほとんど重なるものは省略されます)。
下部はリード分布で、ピーク領域(drompa_peakcallで抽出されるものと同一です)は赤、それ以外の領域は緑で表示されます(この図だとほとんど赤になっていますが)。

1段に表示されるゲノム幅(x軸)は-ls、y軸のスケールは-scale_tag、1ぺージ当たりに表示される段数は-lppで指定することができます。
以下のようにパラメータを変えてみましょう。

$ drompa_draw PC_SHARP -i parse2wigdir/H3K4me3,parse2wigdir/Input,H3K4me3 \
-p K562_H3K4me3_2 -gene refFlat.txt -gt genometable.txt -chr 1 \
-ls 100 -scale_tag 200 -lpp 3

出力は以下のようになります。 f:id:rnakato:20171212184611p:plain

複数サンプルの可視化

4種類のヒストン修飾を並べて可視化してみましょう。 それぞれのサンプル(ChIP・Inputペア)を-iで指定します。コマンドは以下のようになります。

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

-iで指定した順番に並べて可視化されます。 f:id:rnakato:20171212185857p:plain

H3K27me3, H3K36me3はbroad peakであり、H3K4me3に比べるとリードの濃縮が弱いので、縦軸が同じスケールだとわかりづらいですね。 そこで、個別にスケールを変えてみましょう。

$ 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 の項目で最後に指定した数字(200,10,10)が、それぞれのサンプルのy軸のスケールになります。ここで指定したスケールは-scale_tagよりも優先されます。 何故カンマが3つ続くのかは、、次回以降に解説しますので、ここでは置いておきます。

f:id:rnakato:20171212185915p:plain

y軸のスケールが変わり、H3K27me3, H3K36me3の濃縮が見えるようになりました。 H3K36me3はtranscribed stateのマーカー、H3K27me3はsuppressive stateのマーカーで、ゲノム上で排他的に修飾が入ることが知られていますが、それがこの図からもわかります。

以上がDROMPA3での最も基本的な可視化です。これからしばらくは可視化機能について解説していこうと思います。

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

DROMPA3解説その3はピーク抽出(peak calling)です。ピーク抽出とは、ゲノムからreadが有意に濃縮している箇所を網羅的に同定する作業です。

インストール

DROMPA3でのピーク抽出は、drompa_peakcall を使います。 DROMPA3のインストール方法についてはこの記事を参照してください。
単にdrompa_peakcallとタイプすると、以下のヘルプが表示されます。

 $ drompa_peakcall
Usage: drompa_peakcall [--version] <command>
Command: PC_SHARP    peak-calling (for sharp mode)
         PC_BROAD    peak-calling (for broad mode)
         PC_ENRICH   peak-calling (enrichment ratio)

Genome table作成

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

genome tableを作成する - Palmsonntagmorgen

parse2wig

前回のエントリでparse2wigを用いて作成したファイルを使います。 とりあえず下記のコマンドでデータは生成できます。

 $ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeUwTfbs/wgEncodeUwTfbsHelas3CtcfStdAlnRep1.bam
 $ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeUwTfbs/wgEncodeUwTfbsHelas3InputStdAlnRep1.bam
 $ parse2wig -i wgEncodeUwTfbsHelas3CtcfStdAlnRep1.bam -o Ctcf -gt genometable.txt -f BAM
 $ parse2wig -i wgEncodeUwTfbsHelas3InputStdAlnRep1.bam -o Input -gt genometable.txt -f BAM

ピーク抽出

drompa_peakcall は以下の3種類の内部コマンドを持ちます。

  • PC_SHARP: 転写因子などの鋭いピーク(sharp peak)用
  • PC_BROAD: 一部のヒストン修飾など広範囲に分布するピーク(broad peak)用
  • PC_ENRICH: 有意差検定を用いずにシンプルにChIP/Input比のみを考慮するコマンド

ちなみにこの3つのコマンドはデフォルトパラメータが異なるだけで、基本的は一緒だと思ってもらって構いません。ここではPC_SHARPを使いましょう。

下記のようにコマンド含めてタイプすると、オプション一覧が表示されます。

 $ drompa_peakcall PC_SHARP
please specify output prefix.
Usage: drompa_peakcall COMMAND [options] -p <output> -gt <genometable> -i <ChIP>,<Input>
       <ChIP> ChIP sample
       <Input>    Input (control) sample
       <output>: Name of output files
       <genometable>: Tab-delimited file describing the name and length of each chromosome

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

前回生成したファイルを指定して、CTCFのピークを抽出してみましょう。

 $ drompa_peakcall PC_SHARP -i parse2wigdir/Ctcf,parse2wigdir/Input -gt genometable.txt -p CTCF

-iは入力とするデータで、ChIPファイルと対応するInputファイルをカンマで区切って指定します。-gtはGenome table、-pは出力ファイル名です。 このコマンドを実行すると、以下のようなログが標準出力に現れます。

======== PC_SHARP ========
drompa_peakcall version 3.3.0
   output prefix: CTCF
   genometable: genometable.txt
   ChIP: parse2wigdir/Ctcf
   Input: parse2wigdir/Input
   Input format: BINARY
   binsize: 100
   smoothing width: 500 bp
   Peak intensity threshold: 0.00
   ChIP/Input normalization: TOTALREAD
   Enrichment threshold: 0.00
   p-value threshold (internal, -log10): 4.00e+00
   p-value threshold (enrichment, -log10): 4.00e+00
   q-value threshold: 1.00e-03
======================================
sample1: genome: ChIP read=4757329, Input read=14682144, ChIP/Input = 0.32

Peak-calling: chr1..chr2..chr3..chr4..chr5..chr6..chr7..chr8..chr9..chr10..chr11..chr12..chr13..chr14..chr15..chr16..chr17..chr18..chr19..chr20..chr21..chr22..chrX..done.
output peaklist...done.

コマンド実行終了後、カレントディレクトリにCTCF.bed と CTCF.xls という2種類のピークリストが生成されます。この2つのファイルの内容は同一ですが、CTCF.bedはBED形式、CTCF.xlsはより詳細な情報を含むcsv形式という違いがあります。

ピーク抽出の閾値

DROMPA3は以下の五種類のピーク抽出閾値パラメータを持ちます。

  threshold:
       -pthre_internal <double>: p-value threshold for ChIP internal (Poisson, default < 1.0e-04)
       -pthre_enrich   <double>: p-value threshold for ChIP/Input enrichment (binomial, default < 1.0e-04)
       -ethre <double>: IP/Input enrichment threshold (default: 2 for PC_ENRICH, 0 for the others)
       -ipm <double>: read intensity threshold of peak summit (default: 0)
       -qthre <double>: q-value threshold (Bonferroni-Hochberg, default < 0.001)

このうち、-ethre (ChIP/Input比の閾値)、-ipm(ピークの高さの閾値)はPC_SHARPモードではデフォルトでは用いないので、残り3つのパラメータが閾値となります。

有意差検定

DROMPA3は内部で以下の2種類の有意差検定を行います。この2段階検定の方法はPeakSeqのアルゴリズムに基づきます*1

  • ChIP internal: ChIPサンプルにおいて、周辺に対してそのウィンドウにreadが有意に濃縮しているか。(負の二項検定)
  • ChIP/Input enrichment: そのウィンドウにおいで、Inputに対してChIPサンプルのreadが有意に濃縮しているか(二項検定)

図で表すと以下のようになります。

f:id:rnakato:20171211180851p:plain

ChIP internalの検定ではピーク様の濃縮領域を候補領域として全て抽出し、そのうちChIP/Input比が有意であるもの、すなわちInputでは濃縮していない領域のみを抽出します。通常、ChIP/Input enrichmentの閾値の値がピーク数にクリティカルに効いてきますので、ピーク抽出の閾値を変えてみたい場合は、-pthre_enrich オプションの値を色々と変えてみると良いでしょう。
なお、drompa_peakcall はInputサンプル無しでもピーク抽出可能ですが、その場合はChIP/Input enrichmentの検定は省略され、ChIP internalで有意とされた領域全てをピークとして出力します。

このように抽出されたピークリストに対し、さらにBenjamini-Hochberg法による多重検定の補正*2を行います(-qthre)。なお経験的には、ピーク数が数千以上の場合はこの閾値はほとんど関係ありません。。

オプション

  • -binsize: デフォルトでは100bpですが、異なるサイズにすることができます。その場合、parse2wigでそのサイズのデータをあらかじめ生成している必要があります。
  • -sm: スムージングをかける幅です。デフォルトでは500bpです。スムージングをかけると細かいノイズの影響を低減することができますが、ピークの幅が広くなるため、モチーフ解析などをしたい場合は -sm 0としてスムージングをオフにすると良いでしょう。
  • -norm: デフォルトでは、有意差検定を行う際にChIP/Input比が1になるようにread数の正規化を行います。-norm 0とすると正規化をオフにします。
  • -includeYM: デフォルトではY染色体ミトコンドリアのピークは無視されますが、このオプションを追加すると出力されるようになります。

他にもいくつかオプションがあるのですが、それはまた別エントリで解説したいと思います。毎回エントリが長くなり過ぎですね…もうすこしコンパクトにしたいのですが。
とりあえず可視化まではこのまま進もうと思います。

*1:Rozowsky J, Euskirchen G, Auerbach RK, et al. PeakSeq enables systematic scoring of ChIP-seq experiments relative to controls. Nat Biotechnol 2009;27:66–75.

*2:Benjamini-Hochberg法についてはこちらのエントリが参考になります:Benjamini-Hochberg法の原理 - ほくそ笑む