PowerPoint プレゼンテーション

Size: px

Start display at page:

Download "PowerPoint プレゼンテーション"

あゆみむらかわ
9 years ago
Views:

1 平成 28 年度 NGS ハンズオン講習会 ChIP-seq 2016 年 7 月 28 日

4 ChIP-seq とは ChIP(Chromatin Immuno Precipitation) + NGS sequencing クロマチン免疫沈降により濃縮したゲノム領域をシーケンスする手法主な解析対象タンパクとDNAの相互作用ヒストン修飾 Licensed under CC-BY 4.0 Togo picture gallery by DBCLS Amelieff Corporation All Rights Reserved 4

6 ChIP-seq とは input と IP ChIP-seq では免疫沈降のバックグラウンドノイズを削減するためコントロールを使用することが多い免疫沈降 (IP) を行っていないサンプルをコントロールとして使用し検出したピークを抗体に非特異的なものとして取り除くために用いる一般にこのコントロールを input と呼ぶ IP input Licensed under CC-BY 4.0 Togo picture gallery by DBCLS Amelieff Corporation All Rights Reserved 6

8 ChIP-seq 解析の流れ代表的なソフトクオリティコントロール Trimmomatic, fastqc, FASTX_Toolkit... マッピング Bowtie, Bowtie2, bwa... ピーク検出 MACS, MACS2, SICER... ピークアノテーション SnpEff, ChIPpeakAnno... モチーフ探索 rgadem... ChIP-seq 解析の一般的な流れであり全ての ChIP-seq で同一の解析を行うわけではない研究の目的やデータに合わせて最適な解析を設計 Amelieff Corporation All Rights Reserved 8

11 ピーク検出ピーク検出ソフトは IP で濃縮した領域のリードの頂点を検出する本日使用するソフト MACS2( デファクトスタンダード ) 被引用数 2750 件 (2016 年 7 月 20 日時点 ) この他に使用されるソフト SICER 被引用数 425 件 (2016 年 7 月 20 日時点 ) MACS Amelieff Corporation All Rights Reserved 11

公開データの取得今回の解析に必要なデータリファレンスゲノム ( 実行済み ) http://support.illumina.

15 公開データの取得酵母のリファレンスゲノムデータの取得方法 $ wget ftp://igenome:[email protected]/saccharomyces_cerevisiae/ncbi/build3.1/saccha romyces_cerevisiae_ncbi_build3.1.tar.gz $ tar zxvf Saccharomyces_cerevisiae_NCBI_build3.1.tar.gz Saccharomyces cerevisiae のリファレンスゲノムをイルミナの Web ページからダウンロードし解凍 ( 実行済み ) $ ls -l /home/ ユーザ名 /Desktop/amelieff/sacCer3/ : -rwxr-xr-x. 1 root root 月 23 11: genome.fa -rwxr-xr-x. 1 root root 月 23 11: genome.fa.fai -rwxr--r--. 1 root root 月 23 11: mask.gtf -rwxr-xr-x. 1 root root 月 23 11: refgene.txt /home/ ユーザ名 /Desktop/amelieff/Scerevisiae/ の解凍したファイル ( 今回使用するもののみ ) を確認 Amelieff Corporation All Rights Reserved 15

16 公開データの取得 fasta ファイルの中身の確認 $ less /home/ ユーザ名 /Desktop/amelieff/Scerevisiae/genome.fa >chri CCACACCACACCCACACACCCACACACCACACCACACACCACACCACACC CACACACACACATCCTAACACTACCCTAACACAGCCCTAATCTAACCCTG GCCAACCTGTCTCTCAACTTACCCTCCATTACCCTGCCTCCACTCGTTAC CCTGTCCCATTCAACCATACCACTCCGAACCACCATCCATCCCTCTACTT ACTACCACTCACCCACCGTTACCCTCCAATTACCCATATCCAACCCACTG : 1 行目 : コンティグ名 2 行目以降 : 実際の配列情報 q で閲覧を終了する Amelieff Corporation All Rights Reserved 16

公開データの取得データの探し方論文中などから得られたアクセッション番号の ERR1231585

公開データの取得データの探し方 SRA Run Selector でデータを確認する 47Runs 合計

公開データの取得データの探し方 ERR1231585(input) と ERR1231597(sample) を選択コントロール (input)

公開データの取得ダウンロード方法 SRA のダウンロードには SRA-Tools を使用する SRA-Tools(http://ncbi.github.

25 公開データの取得ダウンロード方法 SRA-Tools のインストール本日はデータを用意済みのため実施しません $ wget $ tar xf sratoolkit centos_linux64.tar.gz $ ln -s sratoolkit centos_linux64/bin/prefetch /usr/local/bin/ $ ln -s sratoolkit centos_linux64/bin/fastq-dump /usr/local/bin/ 参考 : Amelieff Corporation All Rights Reserved 25

26 公開データの取得ダウンロード方法 SRA-Tools の prefetch コマンドでまとめて SRA をダウンロードダウンロードした Accession List(SRR_Acc_List.txt) を --option-file で指定 $ prefetch --option-file SRR_Acc_List.txt デフォルトで SRA は ~/ncbi/public/sra/ に保存される $ ls ~/ncbi/public/sra/ ERR sra ERR sra Amelieff Corporation All Rights Reserved 26

27 公開データの取得 SRA の変換方法 SRA-Tools の fastq-dump を使用して SRA から FASTQ へ変換する変換データを保存するディレクトリ (data) を作成する ( 実行済み ) $ mkdir data $ cd data --split-files を付けてペアエンドのファイルを分割しながら FASTQ に変換する ( 実行済み ) どこでペアエンドかシングルエンドかを確認するのか ( 次のスライドで解説 ) $ fastq-dump ~/ncbi/public/sra/err sra --split-files $ fastq-dump ~/ncbi/public/sra/err sra --split-files Amelieff Corporation All Rights Reserved 27

30 公開データの取得実習用データの作成 seqtk( を使用し実習用に FASTQ からデータの一部を抜粋する seqtk のインストール ( 今回は実施しません ) $ wget $ tar xf v1.2.tar.gz $ cd seqtk-1.2 $ ln -s ~/src/seqtk-1.2/seqtk /usr/local/bin/ Amelieff Corporation All Rights Reserved 30

31 公開データの取得実習用データの作成 seqtk を使用し実習用に FASTQ からデータの一部を抜粋する seqtk の実行 $ seqtk sample -s 100 ERR _1.fastq > input_1.fastq $ seqtk sample -s 100 ERR _2.fastq > input_2.fastq $ seqtk sample -s 100 ERR _1.fastq > sample_1.fastq $ seqtk sample -s 100 ERR _2.fastq > sample_2.fastq -s 100: シード値を100に指定ペアで同じシード値を使うことでランダムに抽出するリードのペアを保つ事ができる :50 万リード抽出 Amelieff Corporation All Rights Reserved 31

33 公開データの取得解析対象のシーケンスデータの取得方法 ( 実行済み ) ダウンロード SRA FASTQ 変換 $ prefetch --option-file SRR_Acc_List.txt $ fastq-dump ~/ncbi/public/sra/err12315*.sra --split-files 実習用の軽量なデータを作成 ( 実行済み ) $ seqtk sample -s 100 ERR _1.fastq > input_1.fastq $ seqtk sample -s 100 ERR _2.fastq > input_2.fastq $ seqtk sample -s 100 ERR _1.fastq > sample_1.fastq $ seqtk sample -s 100 ERR _2.fastq > sample_2.fastq Amelieff Corporation All Rights Reserved 33

34 公開データの取得解析対象のシーケンスデータの確認 $ cd /home/iu/chipseq $ ls data input_1.fastq.gz sample_1.fastq.gz input_2.fastq.gz sample_2.fastq.gz アクセッション番号との対応 input_1 ERR _1.fastq.gz input_2 ERR _2.fastq.gz sample_1 ERR _1.fastq.gz sample_2 ERR _2.fastq.gz それぞれ 500,000 リードのデータ Amelieff Corporation All Rights Reserved 34

35 クオリティコントロール QC 前の品質確認シーケンスクオリティチェックソフトウェア FastQC の実行 $ mkdir fastqc_before $ fastqc --nogroup -t 2 -o./fastqc_before data/input_1.fastq.gz data/input_2.fastq.gz data/sample_1.fastq.gz data/sample_2.fastq.gz $ ls fastqc_before input_1_fastqc input_2_fastqc.zip sample_2_fastqc input_1_fastqc.zip sample_1_fastqc sample_2_fastqc.zip input_2_fastqc sample_1_fastqc.zip Amelieff Corporation All Rights Reserved 35

36 クオリティコントロール QC 前の品質確認 FastQC の結果確認 (QC 前 ) 解析結果の html ファイルをブラウザ (firefox) で確認 $ firefox fastqc_before/input_1_fastqc/fastqc_report.html fastqc_before/input_2_fastqc/fastqc_report.html fastqc_before/sample_1_fastqc/fastqc_report.html fastqc_before/sample_2_fastqc/fastqc_report.html ブラウザでタブが 4 つ開かれクオリティチェックの解析結果が確認できる Amelieff Corporation All Rights Reserved 36

クオリティコントロール QC 前の品質確認 fastqc summary (QC 前 ) input1 input2

38 クオリティコントロール QC 処理今回のデータに対する処理 (Trimmomatic を用いた一括処理 1) $ mkdir trimmed_data $ java -jar /usr/local/bin/trimmomatic-0.36.jar PE -threads 2 -phred33 data/input_1.fastq.gz data/input_2.fastq.gz trimmed_data/input_1_paired.fastq trimmed_data/input_1_unpaired.fastq trimmed_data/input_2_paired.fastq trimmed_data/input_2_unpaired.fastq LEADING:20 TRAILING:20 SLIDINGWINDOW:4:15 MINLEN:36 sample~ でも同様の処理を実行 Amelieff Corporation All Rights Reserved 38

39 クオリティコントロール QC 処理今回のデータに対する処理 (Trimmomatic を用いた一括処理 2) $ java -jar /usr/local/bin/trimmomatic-0.36.jar PE -threads 2 -phred33 data/sample_1.fastq.gz data/sample_2.fastq.gz trimmed_data/sample_1_paired.fastq trimmed_data/sample_1_unpaired.fastq trimmed_data/sample_2_paired.fastq trimmed_data/sample_2_unpaired.fastq LEADING:20 TRAILING:20 SLIDINGWINDOW:4:15 MINLEN:36 CPU のコア数に余裕があれば threads の数値を大きくすることでより高速に処理することが可能 Amelieff Corporation All Rights Reserved 39

41 クオリティコントロール QC 後の品質確認 FastQC の結果確認 (QC 後 ) $ mkdir fastqc_after $ fastqc --nogroup -t 2 -o fastqc_after trimmed_data/input_1_paired.fastq trimmed_data/input_2_paired.fastq trimmed_data/sample_1_paired.fastq trimmed_data/sample_2_paired.fastq $ firefox fastqc_after/input_1_paired_fastqc/fastqc_report.html fastqc_after/input_2_paired_fastqc/fastqc_report.html fastqc_after/sample_1_paired_fastqc/fastqc_report.html fastqc_after/sample_2_paired_fastqc/fastqc_report.html Amelieff Corporation All Rights Reserved 41

クオリティコントロール QC 前の品質確認 input_1 input2 sample1 sample2

44 マッピング Bowtie2 によるマッピング (input ファイル ) $ mkdir mapping $ bowtie2 -p 2 -x /home/iu/genome/saccer3/bowtie2index/genome -1 trimmed_data/input_1_paired.fastq -2 trimmed_data/input_2_paired.fastq samtools view -Sb - > mapping/input.bam $ samtools sort mapping/input.bam -o mapping/input.sorted.bam bowtie2 のオプション -p : 使用するスレッド数 -x : bowtie2 で作成したゲノムファイルインデックス -1,-2: 入力 fastq ファイル Samtools のオプション view: SAM もしくは BAM の中身を表示 -Sb: SAM から BAM へ変換 Amelieff Corporation All Rights Reserved 44

45 マッピング Bowtie2 によるマッピング (sample ファイル ) $ bowtie2 -p 2 -x /home/iu/genome/saccer3/bowtie2index/genome -1 trimmed_data/sample_1_paired.fastq -2 trimmed_data/sample_2_paired.fastq samtools view -Sb - > mapping/sample.bam $ samtools sort mapping/sample.bam -o mapping/sample.sorted.bam Amelieff Corporation All Rights Reserved 45

46 ピーク検出 MACS2 によるピーク検出 $ macs2 callpeak -t mapping/sample.sorted.bam -c mapping/input.sorted.bam --outdir macs2_res -f BAMPE -n handson2016 -B -q g 1.2e+7 -t ターゲットサンプル (IP) のファイル -c -t に対するコントロール (input) サンプルのファイル --outdir 結果を出力するディレクトリ -f -t で指定したファイルのファイル形式 BAM SAM BED 他様々なフォーマットが指定可能 BAMPE は paired-end read をマッピングした bam ファイル ( コマンドの説明は次スライドに続きます ) Amelieff Corporation All Rights Reserved 46

47 ピーク検出 MACS2 によるピーク検出 $ macs2 callpeak -t mapping/sample.sorted.bam -c mapping/input.sorted.bam --outdir macs2_res -f BAMPE -n handson2016 -B -q g 1.2e+7 -n 出力ファイルの接頭文字 -B フラグメントの pileup control lambda 値などを BedGraph 形式で保存 -q -g peakcall するピークの閾値 (Benjamini-Hochberg による FDR の q 値 ) デフォルト 0.01 反復領域を除いたゲノムサイズ一部のモデル生物では数字ではなくヒト :hs マウス :mm などの省略が可能 Amelieff Corporation All Rights Reserved 47

48 ピーク検出 MACS2によるピーク検出 $ ls macs2_res handson2016_control_lambda.bdg handson2016_summits.bed handson2016_peaks.narrowpeak handson2016_treat_pileup.bdg handson2016_peaks.xls 各出力ファイルの解説は NGS Surfer s Wikiが参考になるこの後のモチーフ探索にはピークの領域情報が記載された handson2016_peaks.narrowpeak を用いる Amelieff Corporation All Rights Reserved 48

49 ピーク検出先頭の 5 行を確認 $ head -5 handson2016_peaks.narrowpeak chri handson2016_peak_ chrii handson2016_peak_ chriv handson2016_peak_ chriv handson2016_peak_ chriv handson2016_peak_ Amelieff Corporation All Rights Reserved 49

50 ピーク検出 handson2016_peaks.narrowpeak の項目解説列 1: 染色体番号 chri 2: ピーク開始位置 : ピーク終了位置 : ピークの名前 handson2016_peak_1 5: ピークのスコア 46 6: ストランド. 7:fold-change :-log10pvalue :-log10qvalue : ピーク開始位置から頂点までの距離 252 例 Amelieff Corporation All Rights Reserved 50

51 可視化 IGV でピークを確認する検出したピークを IGV で可視化する BAMファイルのインデックスを作成 $ cd mapping $ samtools index input.sorted.bam $ samtools index sample.sorted.bam IGV で下記のファイルを表示 1. handson2016_peaks.narrowpeak 2. input.sorted.bam 3. sample.sorted.bam Amelieff Corporation All Rights Reserved 51

可視化 IGV でピークを確認するスコアの高いピークを確認 $ cd macs2_res $ cat handson2016_summits.

可視化 IGV でピークを確認する検索したピークの範囲が示されている handson2016_peaks.

可視化 IGV でピークを確認するそれぞれのトラックの名前を右クリックし Squished

可視化 IGV でピークを確認する今度は拡大してペアのリードを確認 Expanded

60 アノテーション handson2016_summits.bed に対して snpeff によるアノテーションを実施アノテーション作業用のディレクトリを作成しアノテーション前のファイルを確認 $ mkdir annotation $ cd annotation $ cat../macs2_res/handson2016_summits.bed Amelieff Corporation All Rights Reserved 60

61 アノテーション handson2016_summits.bed に対して snpeff によるアノテーションを実施する $ java -jar /usr/local/bin/snpeff.jar eff -csvstats stats.txt -c /usr/local/bin/snpeff.config -i bed -o bedann R /macs2_res/handson2016_summits.bed > handson2016_summits.annotated.bed eff -csvstats 入力ファイルにアノテーションを行う csv 形式のサマリーファイルを作成する -c snpeff の設定ファイルを指定 -i 入力ファイルのフォーマット -o 出力ファイルのフォーマット ( コマンドの説明は次スライドに続きます ) Amelieff Corporation All Rights Reserved 61

62 アノテーション handson2016_summits.bed に対して snpeff によるアノテーションを実施する $ mkdir annotation $ cd annotation $ java -jar /usr/local/bin/snpeff.jar eff -csvstats stats.txt -c /usr/local/bin/snpeff.config -i bed -o bedann R /macs2_res/handson2016_summits.bed > handson2016_summits.annotated.bed R /macs2_res/handson2016_ summits.bed アノテーションに使用するゲノムバージョン入力ファイル Amelieff Corporation All Rights Reserved 62

63 アノテーション snpeff を用いたアノテーション方法 $ less handson2016_summits.annotated.bed : # Chromo Start End Variant;Annotation Score I I:114304;EXON:ATS1 I I:114304;GENE:YAL019W-A I I:114304;UPSTREAM:FUN30 I I:114304;DOWNSTREAM:LDS1 : 検出されたピークの summit について遺伝子名とその遺伝子に対してエクソン上流下流などの情報が付与される Amelieff Corporation All Rights Reserved 63

64 モチーフ検索 R Bioconductor package rgadem を用いた de novo モチーフ検索 1 $ mkdir../motif $ cd../motif $ R R version ( ) -- "Full of Ingredients" Copyright (C) 2015 The R Foundation for Statistical Computing Platform: x86_64-pc-linux-gnu (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. : Amelieff Corporation All Rights Reserved 64

65 モチーフ検索 R Bioconductor package rgadem を用いた de novo モチーフ検索 2 > library(rgadem) > library("bsgenome.scerevisiae.ucsc.saccer3") > BED <- read.table("../macs2_res/handson2016_peaks.narrowpeak", header=false, sep=" t") > BED <- data.frame(chr=as.factor(bed[,1]), start=as.numeric(bed[,2]), end=as.numeric(bed[,3])) MACS2 から出力された BED ファイルをデータフレームとして読み込む再び handson2016_peaks.narrowpeak を使用 Amelieff Corporation All Rights Reserved 65

66 モチーフ検索 R Bioconductor package rgadem を用いた de novo モチーフ検索 3 > rgbed <- IRanges(start = BED[, 2], end = BED[, 3]) > Sequences <- RangedData(rgBED, space = BED[, 1]) > gadem <- GADEM(Sequences, verbose = 1, genome = Scerevisiae) > pdf("motif.pdf") > plot(gadem) > dev.off() > q() ピーク領域に頻出するモチーフを取得し PDF にプロット Amelieff Corporation All Rights Reserved 66

68 まとめ ChIP-seq 解析の流れクオリティコントロール Trimmomatic, fastqc マッピング Bowtie2 ピーク検出 MACS2 ピークアノテーション SnpEff モチーフ探索 rgadem ChIP-seq 解析の一般的な流れであり全ての ChIP-seq で同一の解析を行うわけではない研究の目的やデータに合わせて最適な解析を設計 Amelieff Corporation All Rights Reserved 68

PowerPoint プレゼンテーション

PowerPoint プレゼンテーション平成 28 年度 NGS ハンズオン講習会 Reseq 解析 2016 年 7 月 26 日本講義にあたって代表的な解析の流れを紹介します論文でよく使用されているツールを使用しますコマンドを沢山実行しますスペルミスが心配な方はコマンド例がありますのでコピーして実行してください実行が遅れてもあせらずに課題や休憩の間に追い付いてください Amelieff Corporation All