大体の解析の流れはこんな感じ。
- 解析に必要なツールのインストール、コントロールを含む.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
- Gene ontology
- BSgenome
- ピーク領域の配列を得る
- rGADEM
- de novo motif search
- MotIV
- 既存データとの類似性チェック
１については既に終わっている。
２から始める。（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があるのでブラウザーにドラッグアンドドロップする。
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.fastqfastqは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したらダウンロードできる。
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とオプションをつける。（結構時間がかかる数十分）
0 件のコメント:
コメントを投稿