大体の解析の流れはこんな感じ。
- 解析に必要なツールのインストール、コントロールを含む.fastqまたは.SRAファイル (SRA Toolkitでfastqへ変換する必要有)をダウンロード
- fastqファイルの中身をFastQCで確認して、FASTX-Toolkitのfastq_quality_trimmer, fastq_quality_filterでクリーニング
- Bowtie2でマッピング(リファレンスゲノムをダウンロードしておく必要有）
- Samtoolsでバイナリー化
- 必要が有ればbedtoolsのintersectBedを使ってリピート配列の除去
- MACSでピークコール
- R+Bioconductorを用いて高次解析
- ChIPpeakAnno
- アノテーション付加（近傍TSSまでの距離・EnsEMBL Gene ID・上流下流）
- org.Mm.eg.db
- BSgenome
- rGADEM
- MotIV
Rを使った高次解析については
二階堂さんのHPが参考になる。とりあえず、1-6に関して実際にやった流れを記す。長いのと、Rの使い勝手がすこしちがうので7に関しては別に記す。
１については既に終わっている。
２から始める。（workディレクトリにfastqファイルがある前提でworkで作業を行う。）
とりあえずfastqc
fastqc -t 8 Sox2_ChIP-seq_Oct3_4_KODay0.fastq -o ./
しばらく待つと、Sox2_ChIP-seq_Oct3_4_KODay0.fastqcというディレクトリとそのzipファイルが出来る。それをディレクトリごとローカルに転送する。ローカル側から
scp -r kosugi@サーバーアカウント:/home/kosugi/local/work/Sox2_ChIP-seq_Oct3_4_KODay0.fastqc.zip /Users/TK/fastqc/
すると指定した場所にディレクトリごと転送される。その中にfastqc_report.htmlがあるのでブラウザーにドラッグアンドドロップする。
こんな感じになる。3'末のクオリティーが低いものが多い事がわかる。FASTQ形式についてwikiで調べればこの数値についてある程度解説が見られる。とりあえず今回は3'末端からクオリティー値を20未満のものをトリミングして、30bp未満のreadを破棄。さらに、80%以上がクオリティー値が20以上ものを抽出する事にした。オプションの詳細はFASTX-Toolkitのusageを確認する。パイプ（｜）でつなげてそのまま実行する。
fastq_quality_trimmer -t 20 -l 30 -Q 33 -i Sox2_ChIP-seq_Oct3_4_KODay0.fastq | fastq_quality_filter -q 20 -p 80 -Q 33 -o ./Sox2_ChIP-seq_Oct3_4_KODay0.clean.fastq
このfastqファイルを先ほどと同様にfastqcしてfastqc_report.htmlをみる。
クリーニングされた事が見た目でわかる。
それぞれのfastqファイルをwc -lすると行数がわかる
wc -l Sox2_ChIP-seq_Oct3_4_KODay0.fastq
62931980 Sox2_ChIP-seq_Oct3_4_KODay0.fastq
wc -l Sox2_ChIP-seq_Oct3_4_KODay0.clean.fastq
59910120 Sox2_ChIP-seq_Oct3_4_KODay0.clean.fastq
fastqは4行で1readなので大体0.8M read程除去されたことがわかる。
3のBowite2によるマッピング
とりあえずコマンド。１０分程かかる。read数やマシーン性能に依存する。
bowtie2 -t -p 8 -x ../bowtie2-2.2.1/indexes/mm10 -U Sox2_ChIP-seq_Oct3_4_KODay0.clean.fastq -S SoxDay0.sam
オプションの意味は
時間：-t
コア数：-p
インデックスのパス：-x
fastqのパス：-U
出力名：-S
4のSamtoolsによるバイナリ化
とりあえずコマンド。
samtools view -Sb SoxDay0.sam > SoxDay0.bam
[samopen] SAM header is present: 66 sequences.等と出て、数分待つ。
オプションの意味は
-S:インプットがsamファイル
-b:アウトプットがbamファイル
bowtieで吐き出されるのは人間が理解できるsamファイルであり、これ以降のステップではバイナリー化しないと認識してくれないらしくこのステップが必要になる。つぎにソートを行う。igvというツールを使う時にindexを作る必要が有り、そのためには順序よく並んでないと怒られるらしい。
samtools sort SoxDay0.bam SoxDay0.sorted
すると、SoxDay0.sorted.bamというのがつくられる。次にindexを作る。
samtools index SoxDay0.sorted.bam
数分待つとSoxDay0.sorted.bam.baiファイルが出来る。今回はigvでマップされたものを確認しない。次にマップ率を出す。flagstatを使う。
samtools flagstat SoxDay0.sorted.bam > SoxDay0.sorted.bam.summary.txt
less SoxDay0.sorted.bam.summary.txt
14977530 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
14529891 + 0 mapped (97.01%:-nan%)
0 + 0 paired in sequencing
0 + 0 read1
0 + 0 read2
0 + 0 properly paired (-nan%:-nan%)
0 + 0 with itself and mate mapped
0 + 0 singletons (-nan%:-nan%)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
97%がマップされている。
5のintersectBedを使ってリピート配列の除去
必要が有ればrepeat region にマップされたreadsを除去する。これには準備が必要で
UCSC Table Browserからデータをダウンロードしておく。以下のように設定して、get outputしたらダウンロードできる。
local/shareに転送した。これで準備完了
scp /Users/TK/Downloads/repeats_mm10.bed kosugi@サーバー名:/home/kosugi/local/share
以下がrepeat regionを除去するコマンド。
intersectBed -abam SoxDay0.sorted.bam -b ../share/repeats_mm10.bed -v > SoxDay0.rmRepeat.bam
これで除去できた。
６のMACSによるピークコール
コントロールのサンプルについても同様の作業を行って、bamファイルが2つある状況にする。準備ができたら以下のコマンドを叩く。数分待つ。
macs14 -t SoxDay0.sorted.bam -c SoxKODay2.sorted.bam -f bam -g mm -n SoxDay0
オプションの意味は
-t:sampleファイル
-c:controlファイル
-f:ファイル形式
-g:ゲノムサイズ
-n:出力名
終わると、色んなファイルが出力される。
SoxDay0_negative_peaks.xls
SoxDay0_summits.bed
SoxDay0_peaks.bed
SoxDay0_peaks.xls
SoxDay0_model.r
例えば、SoxDay0_peaks.bedの中身はこんな感じ。
head SoxDay0_peaks.bed
chr1 3062864 3063136 MACS_peak_1 81.45
chr1 3482902 3483152 MACS_peak_2 65.14
chr1 4150841 4151273 MACS_peak_3 222.85
chr1 4660129 4660269 MACS_peak_4 53.63
chr1 4802472 4802926 MACS_peak_5 193.69
chr1 4855364 4855775 MACS_peak_6 270.83
chr1 4972272 4972486 MACS_peak_7 71.27
chr1 4972497 4972701 MACS_peak_8 84.01
chr1 6448821 6449289 MACS_peak_9 126.84
chr1 6458167 6458317 MACS_peak_10 62.50
以下のように、ローカル側からxlsファイルを転送すればエクセルで見られる。3000ヶ所くらいピークが見つかった。
scp kosugi@サーバー名:/home/kosugi/local/work/SoxDay0_peaks.xls /Users/TK/MACS/
wiggleファイルが欲しい場合はMACSを実行する時に -wとオプションをつける。（結構時間がかかる数十分）
