Palmsonntagmorgen

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

DROMPA3: その2 parse2wig

ChIP-seq DROMPA NGS

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で、約４割がフィルタされています。かなりフィルタされていますね。 なお、値がカッコで囲われているのは、マップリード総数（8,001,988）が指定された数（１０００万）より少ないためです。
その下の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塩基単位。

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