Palmsonntagmorgen

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

記事一覧

解析環境構築

環境変数PATHの通し方

pyenvでPython環境を構築する - Palmsonntagmorgen

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

データ生成

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

2bit genome を作成する - Palmsonntagmorgen

genome tableを作成する - Palmsonntagmorgen

 

fastqデータ取得

SRAからfastqを取得する - Palmsonntagmorgen

ENA,DDBJからfastqを取得する - Palmsonntagmorgen

 

ゲノムマッピング

Readをゲノムにマッピング (その1) - Palmsonntagmorgen

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

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

 

マップデータ操作

SAMtoolsとリダイレクト

SAMtoolsワンライナー覚書

 

ChIP-seq解析: DROMPA

DROMPA3: その1 インストール

DROMPA3: その2 parse2wig

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

 

ChIP-seq解析の注意点

Library complexity (PCR bias)とは何か

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法の原理 - ほくそ笑む

Library complexity (PCR bias)とは何か

前回の DROMPA3: その2 parse2wig - Palmsonntagmorgen で登場した評価指標である Library complexity (PCR bias) の解説です。

PCR biasとは

クロマチン免疫沈降法(ChIP)で得られたリードをゲノムにマップすると、以下のようなマップリード分布が得られるでしょう。 f:id:rnakato:20171126162739p:plain

この図ではゲノム上(下部の黒い線)にリードがマップされている様子を示しています。ゲノム上にはリピート領域やCopy number variation (CNV)などがあるためリード分布は完全に一様にはなりませんが、それでもある程度ランダムにゲノム上に張り付きます。 リードがたくさんマップされている(濃縮している、enrichしている、とも言います)領域は、目的のタンパクが結合する領域かもしれません。

それでは、リードが以下のような分布をしていたらどうでしょうか。

f:id:rnakato:20171126162740p:plain

リード数は最初の例と同じです。このような分布をしていた場合、リードが張り付いている領域をタンパク結合領域だと考えてもよいのでしょうか。
結論から言えば、この例のような場合、リードがたくさんマップされていても濃縮領域であるとは考えません。

ヒトのChIP-seqサンプルのリード数は通常多くとも3000万リード程度ですが、ヒトゲノムは約30億塩基対ありますので、マップリードが一様に分布しているとすると、100塩基対あたり1リードがマップされる計算になります。順鎖、逆鎖に5割ずつマップされていることを考慮すれば、鎖あたりでは200塩基対あたり1リードになります。
その程度の確率の中で、上図のように全く同じポジションに多量のリードが張り付くというのは偶然の産物とは考えられませんので、これらはPCR増幅の産物であると考えられます。

このような偏りのある(PCR biasの大きい)データをそのまま解析に用いると偽陽性(false positives)の原因になりますので、このような同一箇所にマップされたリードは1つを残して全て以下のように除去されます。

f:id:rnakato:20171126182801p:plain

PCR biasをフィルタした結果、この例ではこの領域に4リードしか残りませんでした。PCR biasとしてフィルタされるリードのことを「冗長なリード(redundant reads)」、フィルタ後に残るリードを「非冗長なリード(nonredundant reads)」と呼びます。このように、見かけのマップリード数が多くても、実際に解析に用いられる非冗長なリード数が少ない場合、解析に必要なリード数を満たさない可能性がありますので注意が必要です。

なお、このPCR biasのフィルタリングはChIP-seqだけでなく、オープンクロマチン解析やHi-C解析などでも利用されています。RNA-seqではリードが同一のポジションにマップされることはあり得るため、通常用いられません。

Library complexityとは

読まれたサンプル内に存在するPCR biasの割合を定量的に評価するため、ENCODE projectで"Library complexity"という指標が提案されました*1。Library complexityはnonredundant fraction (NRF)、すなわちマップリード数に対する非冗長なリードの割合によって計測されます。

{ \displaystyle
NRF = \frac{非冗長なリードの数}{マップされたリードの数}
}

しかしこれだけだとまだ問題があります。同じサンプルでもリード数が増えるほど偶然同じ場所にマップされるリードの数が増えるため、同一サンプルでもマップリード数によってNRFの値が変わってしまいます。そこでENCODEでは1000万マップリードあたりの非冗長なリードの数をスコアに用いることを提唱しています。

{ \displaystyle
NRF = \frac{1000万マップリードあたりの非冗長なリードの数}{10,000,000}
}

parse2wigのスタッツファイルで出力されるLibrary complexityはこの値を示します。 ヒトのサンプルでは、この値が0.8以上あることが望ましいとされています。

PCR bias の要因

私の経験上、PCR biasが大きい、すなわちLibrary complexityが低くなる原因の多くは、PCR増幅をかける前の初期DNA量が少ないことによります。
たとえば抗体のS/N比が悪く、タンパク結合DNAがほとんど回収できていない、あるいは細胞数が少なすぎるような場合に、PCR増幅を多くかけることによってシーケンシングのためのDNA量を確保しようとしますが、そのようなサンプルは得てしてLibrary complexityが低くなります。そうすると、解析に使えるリード数は結局足りないということになりがちです。更に、Library complexityの非常に低いサンプルはリードのGC含量に強い偏りがあるなど、他の品質評価基準を満たさないことも多いため、リードを読み足すよりはサンプルを作り直す方が望ましいです。

一方、DNAの断片化に超音波破砕ではなくDNaseを用いた場合、制限酵素サイトによってDNA断片の端が切り揃うのでLibrary complexityが低くなるのでは、という疑問が湧く方もいるかもしれません。私の解析したサンプルの中では、体感ではほとんど違いがないか、あるとしても2~3%程度です。ので、断片化手法の違いはあまり気にする必要はありません。

Library complexityについての考察

私のこれまでの経験では、クオリティが良いとされるサンプルは確かにだいたいLibrary complexityの値が0.8以上ありました。が、0.8を切っているサンプルが直ちに使えないという訳ではありません。0.6程度だったとしても、フィルタ後に残った非冗長なリードの絶対数が十分あれば(最低1000万くらい)、問題なく解析できることもあります。
また、RNA pol2のように非常に抗体のS/N比が高いサンプルでは、リードの濃縮領域に多量にマップされた「本物の」リードがある程度フィルタされてしまうため、Library complexityが低くなります。真のマップリードの一部を除去してしまうことの是非には議論がありますが、ピーク抽出の感度への影響は「ほとんど気にしなくてよいレベル」であると報告されています*2
また、酵母のような小さいゲノムの場合、最初に計算した「塩基あたりにマップされるリードの数」が1以上となります。この場合、1カ所につき1リードのみ許すという基準は厳しすぎるものになりますので、基準を多少緩める必要があります。

まとめ

まとめると、Library complexityのスコアと非冗長なリードの絶対数を確認するのは重要である一方、その値が問題なければ、実際の解析においてPCR biasを除去する必要性自体はそこまでない(かもしれない)ということになります。公開されているChIP-seq解析ツールの大半は(DROMPAも含め)PCR biasをデフォルトで除去しますので、敢えてオフにする必要は特にありませんが、S. cerevisiaeのrDNA 領域のようなhighly repetitiveな領域を観察したい場合には、このフィルタリングをオフにする必要があります。

*1:Landt SG, Marinov GK, Kundaje A, et al. ChIP-seq guidelines and practices of the ENCODE and modENCODE consortia. Genome Res 2012;22:1813-31

*2:Chen Y, Negre N, Li Q, et al. Systematic evaluation of factors influencing ChIP-seq fidelity. Nat Methods 2012;9:609–14

DROMPA3: その2 parse2wig

parse2wigはマッピングファイルを入力とし、Wig形式に変換してくれるツールです。 内部でPCR biasのフィルタ、Total readによる正規化、種々の品質評価も行います。

インストール

DROMPA3のインストール方法についてはこの記事を参照してください。
単にparse2wigとタイプすると、バージョン情報とヘルプが表示されます。

 $ parse2wig
parse2wig version 3.2.7
Usage: parse2wig [option] {-pair} -i <inputfile> -o <output> -gt <genome_table>

       <inputfile>    Mapping file. Multiple files are allowed (separated by ',')
       <output>   Prefix of output files
       <genome_table> Tab-delimited file describing the name and length of each chromosome

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

Genome table作成

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

genome tableを作成する - Palmsonntagmorgen

BAMファイルダウンロード

実際に使ってみた方がわかりやすいと思うので、ここではENCODE projectに登録されているhumanのCTCF ChIP-seqデータ(HeLa S3細胞)を例としてparse2wigを使ってみましょう。
以下のコマンドでBAMファイルをダウンロードします。

 $ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeUwTfbs/wgEncodeUwTfbsHelas3CtcfStdAlnRep1.bam
 $ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeUwTfbs/wgEncodeUwTfbsHelas3InputStdAlnRep1.bam

ここではCTCFサンプルの他に、対応するInputサンプル(ネガティブコントロール)も一緒にダウンロードしています。 初期のChIP-seq解析ではChIPサンプルのみを用いてピーク抽出していましたが、それだとリピート領域やcopy numberの多い領域なども全てピークとして拾ってきてしまうことになるので、 現在はInputのリード分布をコントロールに用いるのが一般的です。

実行

ダウンロードしたBAMファイルを以下のようにしてparse2wigでパースします。

 $ parse2wig -i wgEncodeUwTfbsHelas3CtcfStdAlnRep1.bam -o Ctcf -gt genometable.txt -f BAM

-i で入力ファイルを指定、-oで出力ファイル名を指定します。-gtでGenome tableを指定します。 -fは入力ファイルのフォーマットで、ここではBAMを指定します。

parse2wigを実行すると、最初に以下のような指定オプション一覧が表示されますので、予定通りのオプションになっているか確認することができます。

======================================
parse2wig version 3.2.7

Input file: wgEncodeUwTfbsHelas3CtcfStdAlnRep1.bam
    Format: BAM, SINGLE_END
Output file: parse2wigdir/Ctcf
    Format: BINARY
Genome-table file: genometable.txt
Binsize: 100 bp
Fragment length: 150
PCR bias filtering: ON
Read number for library complexity: 10.0 M
Cross-correlation profile: OFF

Total read normalization: NONE
======================================

Binsizeはリードをカウントするウィンドウサイズです。 マップされたリードはFragment lengthで指定された値まで伸長してカウントされます。
PCR bias filteringは、ゲノムの同一部位にマップされた複数のリードをPCR biasとして除去する作業で、 ChIP-seq解析では一般的に行う作業ですが、リピート領域などを特に解析したい場合にはオフにすることもあります。
Library complexityは、入力ファイルに含まれているマップリードの数のうち、PCR bias除去後に残ったリード数の割合を意味します。 この値は入力リード数に依存しますので、1000万マップリードにおけるPCR bias除去後に残ったリード数の割合を出力します。
Cross-correlation profileは抗体のS/N比を示す値で、いずれ別記事で説明しようと思います。
Total read normalizationは、マップリード総数でデータを正規化する作業です。得られるピークの数はマップリード数に依存しますので、サンプル間でピーク数を比較するような場合はこの正規化が必要になります。

同様に、Inputサンプルもparse2wigにかけましょう。

 $ parse2wig -i wgEncodeUwTfbsHelas3InputStdAlnRep1.bam -o Input -gt genometable.txt -f BAM

出力

parse2wigを実行すると、parse2wigdirというディレクトリが自動的に作成され、その中に出力ファイルが生成されます。 ls でparse2wigdirを表示すると以下のようになっているはずです。

 $ ls parse2wigdir/
Ctcf.100.binarray_dist.xls  Input.100.binarray_dist.xls
Ctcf.100.xls                Input.100.xls
Ctcf.readlength_dist.xls    Input.readlength_dist.xls
Ctcf_chr1.100.bin           Input_chr1.100.bin
Ctcf_chr10.100.bin          Input_chr10.100.bin
Ctcf_chr11.100.bin          Input_chr11.100.bin
Ctcf_chr12.100.bin          Input_chr12.100.bin
Ctcf_chr13.100.bin          Input_chr13.100.bin
Ctcf_chr14.100.bin          Input_chr14.100.bin
Ctcf_chr15.100.bin          Input_chr15.100.bin
Ctcf_chr16.100.bin          Input_chr16.100.bin
Ctcf_chr17.100.bin          Input_chr17.100.bin
Ctcf_chr18.100.bin          Input_chr18.100.bin
Ctcf_chr19.100.bin          Input_chr19.100.bin
Ctcf_chr2.100.bin           Input_chr2.100.bin
Ctcf_chr20.100.bin          Input_chr20.100.bin
Ctcf_chr21.100.bin          Input_chr21.100.bin
Ctcf_chr22.100.bin          Input_chr22.100.bin
Ctcf_chr3.100.bin           Input_chr3.100.bin
Ctcf_chr4.100.bin           Input_chr4.100.bin
Ctcf_chr5.100.bin           Input_chr5.100.bin
Ctcf_chr6.100.bin           Input_chr6.100.bin
Ctcf_chr7.100.bin           Input_chr7.100.bin
Ctcf_chr8.100.bin           Input_chr8.100.bin
Ctcf_chr9.100.bin           Input_chr9.100.bin
Ctcf_chrM.100.bin           Input_chrM.100.bin
Ctcf_chrX.100.bin           Input_chrX.100.bin
Ctcf_chrY.100.bin           Input_chrY.100.bin

*.binが生成されたWigファイルで、デフォルトではBinary形式になります。Binary形式の場合、染色体ごとにファイルが出力されます。 *.binarray_dist.xlsは各ウィンドウにマップされたリード数の分布、*.readlength_dist.xlsはリード長の分布を保存しています。
ここで生成されたデータを以降の解析で用います。

スタッツ確認

parse2wigdir/Ctcf.100.xls には入力サンプルのクオリティチェックの結果が記載されています。 エクセルなどで表示してみましょう。ここではlibreoffice --calcを使います。

 $ libreoffice --calc parse2wigdir/Ctcf.100.xls

f:id:rnakato:20171116164411p:plain

スタッツファイルの中身はこのようになっています。 Library complexityは 0.595で、約4割がフィルタされています。かなりフィルタされていますね。 なお、値がカッコで囲われているのは、マップリード総数(8,001,988)が指定された数(1000万)より少ないためです。
その下のPossionとNegative binomialはピークコールをするために推定されたパラメータです。

その下に、全ゲノム及び各染色体のパラメータが並びます。内訳は以下の通りです。

  • length: 配列長
  • mappable base: マップされ得る配列長(ここではマッパビリティファイルを指定していないので、配列長と同じになっています)
  • mappability: lengthに対するmappable baseの割合
  • total reads: マップされた数(forward, reverse, bothの三種類)
  • %genome: 各染色体にマップされたリード数の全リード数に対して占める割合
  • nonredundant reads: PCR bias filtering後に残ったリードの数。PCR bias filteringをオフにするとtotal readsと同じになります
  • redundant reads: PCR bias としてフィルタされたリードの数
  • read depth: ゲノム長に対するリード数の多さ。(リード数 x fragment length / genome length)
  • scaling weight: total read normalization をオンにした場合に、正規化係数が表示される
  • normalized read number: 正規化後のリード数。ここでは正規化がオフなので、nonredundant readsと値が同じ
  • genome coverage: mappable base に対する、リードがマップされた塩基の割合。1塩基単位。

一度に書いてもわかりにくいと思いますので、次回もう少し詳細に説明します。

DROMPA3: その1 インストール

今回からは私が開発したDROMPA3の利用法について解説します。

DROMPAとは

ChIP-seq解析のためのパイプラインツールです。ピーク抽出の他、品質評価、可視化、複数サンプルの比較解析などができます。 複数のサンプルを同時に解析できること、pdf形式でデータを出力できること、ChIP-seq解析の様々なステップをオールインワンで含んでいることが特長です。 以下が元論文です。

DROMPA: easy-to-handle peak calling and visualization software for the computational analysis and validation of ChIP-seq data - Nakato - 2013 - Genes to Cells - Wiley Online Library

この論文はversion 1の頃のもので、現在のversion3は色々な部分が変わっているので、詳細はDROMPA3付属のマニュアルを読んでいただくと良いと思います。
ちなみに名前の由来は、DRaw and Observe Multiple enrichment Profiles and Annotation の頭文字を取ったもの…ということになっていますが、 本当はFC東京のマスコットのドロンパから名前をもらいました。(FC東京サポなので)

以下、インストール方法です。

関連ライブラリのインストール

まずインストールに必要となるライブラリ群をインストールします。詳細はDROMPA3のWebサイトを参照してください。

Ubuntuの場合:

 $ sudo apt install git gcc libgtk2.0-dev libgsl-dev

CentOSの場合:

 $ sudo yum -y install zlib-devel gsl-devel gtk2-devel

ダウンロード・コンパイル

githubからダウンロードし、以下のようにしてコンパイルします。

 $ git clone https://github.com/rnakato/DROMPA3
 $ cd DROMPA3
 $ make

関連ツールのインストール

BAMフォーマットのファイルを入力にする場合、SAMtools が必要になります。 インストール方法は以下の記事を参照してください。

SAMtoolsとリダイレクト - Palmsonntagmorgen

また、複数のpdfを結合するために Coherent PDF というツールを使いますので、こちらもgithubからダウンロードします。

 $ git clone https://github.com/coherentgraphics/cpdf-binaries

PATHを通す

最後にPATHを通します。PATHとは何ぞや?という方は以前の記事を参照してください。
例としてここではホームディレクトリに my_chipseq_exp というディレクトリを作り、その中にDROMPA3とCoherent PDFをダウンロードした場合のPATHを記載します。

 $ export PATH = $PATH:$HOME/my_chipseq_exp/DROMPA3:$HOME/my_chipseq_exp/cpdf-binaries/Linux-Intel-64bit/

これでインストールは完了です。 以下のコマンドをタイプして、DROMPA3のヘルプが表示されれば成功です。

 $ 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
         FRIP        accumulate read counts in bed regions specified
         CI          compare peak-intensity between two samples
         CG          output ChIP-reads in each gene body
         GOVERLOOK   genome-wide overlook of peak positions
         PROFILE     make R script of averaged read density
         HEATMAP     make heatmap of multiple samples
         TR          calculate the travelling ratio (pausing index) for each gene

SAMtoolsワンライナー覚書

順次追加するかも。versionは1.5です。
.sort.bam はソート済BAMを表します。

SAM -> BAM 変換

 $ samtools view -bS sample.sam > sample.bam

BAM -> SAM 変換

 $ samtools view sample.bam > sample.sam

BAMをソート

 $ samtools sort sample.bam > sample.sort.bam

SAM -> ソート済BAM

 $ samtools view -bS sample.sam | samtools sort > sample.sort.bam
 または
 $ samtools sort sample.sam > sample.sort.bam

BAM indexを作成

 $ samtools index sample.sort.bam

複数のBAMファイルをマージ

 $ samtools merge sample.merged.bam sample1.bam sample2.bam

SAM/BAMに含まれるリード数をカウント(unmappedも含む)

 $ samtools view -c sample.bam

SAM/BAMに含まれるリード数をカウント(mapped readのみ)

 $ samtools view -F 0x04 -c sample.bam

マップされなかったリード数をカウント

 $ samtools view -f4 -c sample.bam

染色体別のマップリード数をカウント(index要)

 $ samtools idxstats sample.bam

もうちょっと詳細な情報が知りたい

 $ samtools flagstat sample.bam

さらに詳細な情報が知りたい

 $ samtools stats sample.bam

BAMからmapped readのみを抽出

 $ samtools view -F 0x04 -b sample.bam > sample.mapped.bam

BAMからPCR duplicateを除去

 $ # single-end
 $ samtools rmdup -s sample.sort.bam sample.sort.rmdup.bam  
 $ # paired-end
 $ samtools rmdup sample.sort.bam sample.sort.rmdup.bam    

BAMからPCR duplicate除去済のmapped readを抽出

 $ samtools view -F 0x04 -b sample.sort.bam | samtools rmdup -s - sample.rmdup.sort.bam