2016_RNAseq解析_修正版

Size: px

Start display at page:

Download "2016_RNAseq解析_修正版"

まなながだき
5 years ago
Views:

1 平成 28 年度 NGS ハンズオン講習会 RNA-seq 解析 2016 年 7 27

5 RNA-seq とは n メッセンジャー RNA(mRNA) をキャプチャして次世代シーケンサーでシーケンシングする法 n リファレンスがある物種の場合 : 既知遺伝にマッピングするリファレンスにマッピングして遺伝発現量を定量する n リファレンスがない物種の場合 : アセンブリングして転写物構造を予測しそれに対してマッピングする近いゲノムのリファレンスにマッピングする Amelieff Corporation All Rights Reserved 5

7 RNA-seq 解析の流れ解析ソフト例クオリティコントロール PRINSEQ, cutadapt, FASTX_Toolkit... マッピング HISAT2, TopHat2, STAR, Bowtie2... 発現定量 Cufflinks, StringTie, RSEM... 発現較 Cuffdiff, DESeq, edger... 可視化 CummeRbund... n RNA-seq 解析の般的な流れであり全ての RNA-seq で同の解析をうわけではありません n 研究の的やデータに合わせて最適な解析をいます Amelieff Corporation All Rights Reserved 7

9 公開データの取得酵のリファレンスゲノムデータの取得法 $ wget ftp://igenome:g3nom3s4u@ussdftp.illumina.com/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/iu/genome/saccer3/ : -rwxr-xr-x. 1 iu iu :09 genome.fa : -rwxr-xr-x. 1 iu iu :09 genome.fa.fai : /home/iu/genome/saccer3/ に解凍したファイル ( 今回使するデータのみ ) を置いてあるので確認する Amelieff Corporation All Rights Reserved 9

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

公開データの取得解析対象のシーケンスデータの取得法 1 http://trace.ddbj.nig.ac.jp/dra/index.

公開データの取得解析対象のシーケンスデータの取得法 3 研究概要を確認する click!

15 公開データの取得解析対象のシーケンスデータの取得法 5 データをダウンロードする ( 実済み ) $ mkdir p rnaseq/data 解析をう rnaseq ディレクトリと data ディレクトリを作成する $ cd rnaseq/data/ $ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/sralite/byexp/lites ra/srx/srx104/srx /srr /srr sra SRR SRR SRR についても同様にダウンロード SRR : RPL_yeast_imidazole_ligase SRR : RPL_yeast_imidazole_ligase_rep2 SRR : RPL_yeast_imidazole_noligase SRR : RPL_yeast_imidazole_noligase_rep2 Amelieff Corporation All Rights Reserved 15

16 公開データの取得解析対象のシーケンスデータの取得法 6 SRA データを FASTQ ファイルへ変換する ( 実済み ) $ fastq-dump --split-files SRR sra fastq-dump コマンドは NCBI SRA toolkit をインストールすると利できる $ head SRR _1.fastq > 10K_SRR _1.fastq $ head SRR _2.fastq > 10K_SRR _2.fastq 先頭 1 万リードを抽出する ( 実済み ) SRR SRR SRR についても同様に処理する $ ls data 10K_SRR _1.fastq 10K_SRR _1.fastq 10K_SRR _2.fastq 10K_SRR _2.fastq 10K_SRR _1.fastq 10K_SRR _1.fastq 10K_SRR _2.fastq 10K_SRR _2.fastq Amelieff Corporation All Rights Reserved 16

17 公開データの取得解析対象のシーケンスデータの取得法 7 シーケンスデータを確認する $ less NS500272:29:H2KGHBGXX:1:11101:4753:1025 length=80 NTGGTNCCGAAGCTCCCACTTATTCTACACCCTCTATGTCTCTTCACAATGTCAAACTAGAGTC AAGCTCAACAGGGTCT +SRR NS500272:29:H2KGHBGXX:1:11101:4753:1025 length=80 #AAAA#FFFFFFFFFFFFFFFFFFFFFF<FFFFFFFFFFFFFFFFFFFFFFFFFFFFAFFFFFF FAFFFFF.FFFFFFFF : fastq ファイルの中を表する 1 配列 IDと付加情報 2 : 塩基配列 3 : + 配列 IDと付加情報 4 : クオリティ fastq ファイルは 1 リードあたり 4 で表記される Amelieff Corporation All Rights Reserved 17

18 公開データの取得解析対象のシーケンスデータの取得法 8 リード数を確認する $ wc -l data/10k_srr _1.fastq data/10k_srr _1.fastq 40,000 が抽出されていることを確認する $ wc -l data/10k_srr _2.fastq data/srr _2.fastq 1リードは4 なのでリード数は 40,000 / 4 = 1 万リードである $ wc l data/* ワイルドカード (*) でまとめて確認できる Amelieff Corporation All Rights Reserved 18

19 クオリティコントロール FastQC: シーケンスクオリティチェックソフトウェア $ fastqc -v FastQC v バージョンを確認する ( 最新版は v0.11.5) $ fastqc -h SYNOPSIS FastQC - A high throughput sequence QC analysis tool fastqc seqfile1 seqfile2.. seqfilen fastqc [-o output dir] [--(no)extract] [-f fastq bam sam] [-c contaminant file] seqfile1.. seqfilen :.fastq 以外に.sam や.bam も指定可能複数ファイルの指定も可能である Amelieff Corporation All Rights Reserved 19

20 クオリティコントロール FastQC の実 $ mkdir fastqc_res $ fastqc -o fastqc_res -f fastq --nogroup data/10k_srr _1.fastq data/10k_srr _2.fastq $ ls fastqc_res 10K_SRR _1_fastqc 10K_SRR _1_fastqc.zip 10K_SRR _2_fastqc 10K_SRR _2_fastqc.zip 解析結果の html ファイルをブラウザ (firefox) で確認する $ firefox fastqc_res/10k_srr _1_fastqc/fastqc_report.html fastqc_res/10k_srr _2_fastqc/fastqc_report.html WEB ブラウザ上でクオリティチェックの解析結果が確認できる SRR SRR SRR についても同様に処理する Amelieff Corporation All Rights Reserved 20

22 クオリティコントロール FastQC の結果確認 2 Per Base Sequence Quality 横軸はリード縦軸は quality value を表すリードの位置における全体のクオリティの中央値や平均を確認できる線は中央値線は平均値のボックスは 25% 75% の領域を表す上下に伸びたいバーが 10% 90% の領域を意味する Amelieff Corporation All Rights Reserved 22

クオリティコントロール FastQC の結果確認 3 Per Sequence Quality Scores 縦軸がリード数横軸が

クオリティコントロール FastQC の結果確認 5 Per Base GC Content リードにおける位置での GC 含量を表すいずれかの位置で全体での GC

クオリティコントロール FastQC の結果確認 6 Per Base N Content N はシーケンサーの問題で ATGC いずれの塩基にも決定出来なかった場合に記述される

クオリティコントロール FastQC の結果確認 7 Sequence Length Distribution リードの全体の分布全てのリードの

クオリティコントロール FastQC の結果確認 8 Sequence Duplication Levels リードの重複レベルをている 1 10 はそれぞれ重複のレベルで全体の 20%

クオリティコントロール FastQC の結果確認 9 Overrepresented Sequences 重複している配列とその割合を表す特定の配列が全リードの

31 クオリティコントロール PRINSEQのインストール ( 実済み ) $ wget seq-lite tar.gz # ソースコードをダウンロード $ tar zxvf prinseq-lite tar.gz # 展開 $ cd prinseq-lite # プログラムのっているディレクトリに移動 $ chmod +x prinseq-lite.pl # 実権限を付与 $ ln s /path/to/prinseq-lite /prinseq-lite.pl /usr/local/bin #PATHを通す PRINSEQの使いの確認 $ prinseq-lite.pl -h Amelieff Corporation All Rights Reserved 31

33 マッピング HISAT2 のインストール ( 実済み ) $ wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/downloads/hisat linux_x86_64.zip # ソースコードをダウンロード $ unzip hisat linux_x86_64.zip # 展開 $ ln s /path/to/hisat /hisat2 /usr/local/bin #PATHを通す Amelieff Corporation All Rights Reserved 33

34 マッピング HISAT2 の Index ファイルのダウンロード ( 実済み ) $ wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/sc3.tar.gz $ tar zxvf sc3.tar.gz $ ls /home/iu/genome/saccer3/hisat2index/ genome.1.ht2 genome.3.ht2 genome.5.ht2 genome.7.ht2 genome.2.ht2 genome.4.ht2 genome.6.ht2 genome.8.ht2 HISAT2 の使いの確認 $ hisat2 -h HISAT2 version by Daehwan Kim (infphilo@gmail.com, Usage:hisat2 [options]* -x <ht2-idx> {-1 <m1> -2 <m2> -U <r> --sra-acc <SRA accession number>} [-S <sam>] : Amelieff Corporation All Rights Reserved 34

35 発現定量遺伝の発現量遺伝上にマップされたリード数 n い遺伝ほどマップされるリードは多くなる ( 遺伝間のバイアス ) n サンプル量の多いランほどマップされるリードは多くなる ( ラン間のバイアス ) これらのバイアスを補正してから発現量を較する必要がある n 発現量としてよく使われる指標 RPKM (Reads Per Kilobase per Million mapped reads) FPKM (Fragments Per Kilobase of exon per Million mapped fragments) どちらも発現量をエクソンと全マッピング数で補正した値 Amelieff Corporation All Rights Reserved 35

36 発現定量 Cufflinks: 遺伝発現解析ソフトウェア cufflinks の使いの確認 $ cufflinks cufflinks v2.2.1 linked against Boost version Usage: cufflinks [options] <hits.sam> General Options: -o/--output-dir write all output files to this directory -p/--num-threads number of threads used during analysis --seed value of random number generator seed -G/--GTF quantitate against reference transcript annotations -g/--gtf-guide -M/--mask-file use reference transcript annotation to guide assembly ignore all alignment within transcripts in this file : -g: アセンブルのガイドとして既知の遺伝情報を使することができる Amelieff Corporation All Rights Reserved 36

37 その他の RNA-seq 解析 n リファレンスがある物種の場合 : SNP InDel 検出 (GATK, freebays) 融合遺伝検出 (Chimerascan, TopHat-Fusion) circular RNA 検出 (DCC, CIRI) シングルセルRNA 解析 (Cell Ranger (10x Genomics)) n リファレンスがない物種の場合 : De novo アセンブリ (Trinity) Amelieff Corporation All Rights Reserved 37

41 クオリティコントロール PRINSEQ によるクオリティコントロール $ mkdir 1_qc $ prinseq-lite.pl -fastq data/10k_srr _1.fastq -fastq2 data/10k_srr _2.fastq -out_good 1_qc/10K_SRR notail -out_bad null -out_format 3 -trim_left 5 -trim_tail_right 5 -trim_qual_right 30 -ns_max_p 20 -min_len 30 -fastq の FASTQ ファイル -fastq2 の FASTQ ファイル ( ペアエンドの場合 ) -out_good フィルターを通過したリードの名前 -out_bad フィルタリングされたリードの名前 (null は出しない ) -out_format 1 (FASTA only), 2 (FASTA and QUAL), 3 (FASTQ), 4 (FASTQ and FASTA), or 5 (FASTQ, FASTA and QUAL) Amelieff Corporation All Rights Reserved 41

42 クオリティコントロール PRINSEQ によるクオリティコントロール $ mkdir 1_qc $ prinseq-lite.pl -fastq data/10k_srr _1.fastq -fastq2 data/10k_srr _2.fastq -out_good 1_qc/10K_SRR notail -out_bad null -out_format 3 -trim_left 5 -trim_tail_right 5 -trim_qual_right 30 -ns_max_p 20 -min_len 30 SRR SRR SRR についても同様に処理する -trim_tail_right -trim_qual_right -ns_max_p 3ʼ 末端のポリテールが 5 以上の末端を除去 3ʼ 末端からクオリティ 30 以下の塩基を除去未知の塩基 (N) が多いリード除去 (20% 以上 ) -min_len 配列が短いリード除去 (30bp 以下 ) PRINSEQ は極めて多機能なソフトウェアでありクオリティチェックからトリミングフィルタリングまで様々なプロセスが可能 Amelieff Corporation All Rights Reserved 42

クオリティコントロール FastQC の実 $ fastqc -o fastqc_res -f fastq --nogroup 1_qc/10K_SRR2048224.notail_1.fastq 1_qc/10K_SRR2048224.notail_2.

43 クオリティコントロール FastQC の実 $ fastqc -o fastqc_res -f fastq --nogroup 1_qc/10K_SRR notail_1.fastq 1_qc/10K_SRR notail_2.fastq 解析結果の html ファイルをブラウザ (firefox) で確認する $ firefox fastqc_res/10k_srr notail_1_fastqc.html クリーニング前クリーニング後 Amelieff Corporation All Rights Reserved 43

48 マッピング $ mkdir 2_mapping $ hisat2 -x /home/iu/genome/saccer3/hisat2index/genome --dta --dta-cufflinks -1 1_qc/10K_SRR notail_1.fastq -2 1_qc/10K_SRR notail_2.fastq -S 2_mapping/10K_SRR sam --dta --dta-cufflinks マッピング結果からアセンブリをう cufflinks のためのアセンブリをう -S SAM ファイルに書き出す名前 SRR SRR SRR についても同様に処理する Amelieff Corporation All Rights Reserved 48

49 マッピング $ mkdir 2_mapping $ hisat2 -x /home/iu/genome/saccer3/hisat2index/genome --dta --dta-cufflinks -1 1_qc/10K_SRR notail_1.fastq -2 1_qc/10K_SRR notail_2.fastq -S 2_mapping/10K_SRR sam 9992 reads; of these: 9992 (100.00%) were paired; of these: 2174 (21.76%) aligned concordantly 0 times 193 (1.93%) aligned concordantly exactly 1 time 7625 (76.31%) aligned concordantly >1 times pairs aligned concordantly 0 times; of these: 3 (0.14%) aligned discordantly 1 time pairs aligned 0 times concordantly or discordantly; of these: 4342 mates make up the pairs; of these: 2925 (67.37%) aligned 0 times 99 (2.28%) aligned exactly 1 time 1318 (30.35%) aligned >1 times 85.36% overall alignment rate Amelieff Corporation All Rights Reserved 49

50 マッピング SAM ファイルを BAM ファイルに変換 $ samtools view -b 2_mapping/10K_SRR sam > 2_mapping/10K_SRR bam $ ls -lh -rw-rw-r-- 1 iu iu 1.7M : K_SRR bam -rw-rw-r-- 1 iu iu 13M : K_SRR sam 13M の SAM ファイルが 1.7M のバイナリファイルに変換される BAM ファイルをソート $ samtools sort 2_mapping/10K_SRR bam -o 2_mapping/10K_SRR sorted.bam $ ls 10K_SRR bam 10K_SRR sam 10K_SRR sorted.bam SRR SRR SRR についても同様に処理する Amelieff Corporation All Rights Reserved 50

51 マッピング結果の可視化 BAM ファイルのインデックスを作成 $ samtools index 2_mapping/10K_SRR sorted.bam $ ls 2_mapping 10K_SRR bam 10K_SRR sam 10K_SRR sorted.bam 10K_SRR sorted.bam.bai : SRR SRR SRR についても同様に処理する Amelieff Corporation All Rights Reserved 51

52 マッピング結果の可視化 Integrative Genomics Viewer (IGV) をいた解析結果の確認 1 $ igv.sh IGV を起動し Genomes タブから Load Genomes from File... を選択 /home/iu/geno me/saccer3 の下にある genome.fa を選択し開く Amelieff Corporation All Rights Reserved 52

マッピング結果の可視化 Integrative Genomics Viewer (IGV) をいた解析結果の確認 3

56 発現定量 cufflinks の実 $ cufflinks -o SRR min-frags-per-transfrag 2 2_mapping/10K_SRR sorted.bam $ ls SRR genes.fpkm_tracking skipped.gtf isoforms.fpkm_tracking transcripts.gtf genes.fpkm_tracking の確認 $ less genes.fpkm_tracking 4 列に Gene ID 10 列に FPKM が記載されている SRR SRR SRR についても同様に処理する Amelieff Corporation All Rights Reserved 56

57 発現定量 transcripts.gtf.txt の作成 $ vim transcripts.gtf.txt 挿モード (i) で以下を記 SRR /transcripts.gtf SRR /transcripts.gtf SRR /transcripts.gtf SRR /transcripts.gtf 挿モードの終了 : エスケープ (ESC) コマンドモード : コロン (:) 保存 : コマンドモードで w + Enter 終了 : コマンドモードで q + Enter 保存せずに終了 : コマンドモードで q! + Enter Amelieff Corporation All Rights Reserved 57

59 発現較 cuffmerge の実 $ cuffmerge -o COMPARE -g /home/iu/genome/saccer3/genes.gtf -s /home/iu/genome/saccer3/genome.fa transcripts.gtf.txt -o/--output-dir 出ディレクトリ -g/--ref-gtf アノテーションのgtfファイル -s/--ref-sequence リファレンスゲノムFASTA ファイル -p/--num-threads スレッド数 ( デフォルト=1) Amelieff Corporation All Rights Reserved 59

60 発現較 cuffdiff の実 $ cuffdiff -o COMPARE -L Group1,Group2 COMPARE/merged.gtf 2_mapping/10K_SRR sorted.bam, 2_mapping/10K_SRR sorted.bam 2_mapping/10K_SRR sorted.bam, 2_mapping/10K_SRR sorted.bam -o/--output-dir アノテーションのgtfファイル -L/--labels グループの指定 ( カンマ区切り ) -p/--num-threads スレッド数 ( デフォルト=1) Amelieff Corporation All Rights Reserved 60

62 可視化 cummerbund の紹介 Cufflinks の結果をいて可視化をう R パッケージ $ R #Cufflinks の実ディレクトリ (COMPARE) で起動する > library(cummerbund) > cuff <- readcufflinks() > cuff CuffSet instance with: 2 samples # サンプル数 6935 genes # 遺伝数 7077 isoforms # 転写産物数 7052 TSS # 転写開始位置数 6643 CDS # コード領域数 6935 promoters # プロモーター数 7052 splicing # スプライシング領域数 6534 relcds # 調節コード領域 Amelieff Corporation All Rights Reserved 62

可視化 cummerbund の紹介 > s <- csscatter(genes(cuff), "Group1 ", "Group2 ", smooth=t)

64 可視化 cummerbund の紹介 > dens <- csdensity(genes(cuff)) > dens > densrep <- csdensity(genes(cuff),replicates=t) > densrep グループごとの log 10 (fpkm) の分布サンプルごとの log 10 (fpkm) の分布 Amelieff Corporation All Rights Reserved 64

最後に本講義でご紹介した流れは解析法の例ですツールの選択に正解はありませんのデータに適したツールを選択しより良い解析順を確

PowerPoint プレゼンテーション

PowerPoint プレゼンテーション V1 次世代シークエンサ実習 II 本講義の内容 Reseq 解析 RNA-seq 解析公開データ取得クオリティコントロールマッピング変異検出公開データ取得クオリティコントロールマッピング発現定量 FPKM を算出します 2 R N A - s e q とはメッセンジャー RNA(mRNA) をキャプチャして次世代シーケンサーでシーケンシングする手法リファレンスがある生物種の場合