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コマンドでのマクロな可視化

ChIP-seq解析の注意点

Library complexity (PCR bias)とは何か

ピークを入力とする操作

BEDtoolsワンライナー覚書

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++19とどんどん進化していっているんだけど、そういう最先端の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 d1.bed d2.bed d3.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の一番の弱点はコマンドが長くなってしまうことなのですが、このようにするとすっきりと記述することができますので、試してみてください。