V1RNA-seq

Copyright © Amelieff Corporation All Rights Reserved.

本講義の内容• Reseq解析• RNA-seq解析

公開データ取得

公開データ取得↓↓

クオリティコントロール

クオリティコントロール↓↓マッピングマッピング

↓↓

変異検出

発現定量
FPKMを算出します。

Copyright © Amelieff Corporation All Rights Reserved.

2RNA-seqとは•メッセンジャーRNA（mRNA）をキャプチャして次世代シーケンサーでシーケンシングする手法•

リファレンスがある生物種の場合：
本日の内容– 既知遺伝子にマッピングする– リファレンスにマッピングして遺伝子構造を同定する•

リファレンスがない生物種の場合：– アセンブリングして転写物構造を予測し、それに対してマッピングする– 近いゲノムのリファレンスにマッピングする

Copyright © Amelieff Corporation. All Rights Reserved.

3RNA-seq解析：パイプライン
データ取得 → クオリティコントロール → マッピング→発現定量

スプライスジャンクションを考慮したマッピングを行います。
Copyright © Amelieff Corporation All Rights Reserved.

4RNA-seq解析：パイプライン
データ取得 → クオリティコントロール → マッピング→発現定量

今日は一部のコマンドを実行します。

ヒートマップ、GO解析、パスウェイ解析 など様々。
Copyright © Amelieff Corporation All Rights Reserved.

5RNA-seq解析でできること• 発現量の定量・比較• 新規転写物・新規スプライシングバリアントの探索

RNA-seqがマイクロアレイと比較して優れている点•••

新規転写物や融合遺伝子が検出可
SNV・small Indelも検出可
プローブの設計を必要としない（非モデル生物にも対応可）

Copyright © Amelieff Corporation. All Rights Reserved.

6RNA-seq解析：データ
データ取得 → クオリティコントロール → マッピング→発現定量• 酵母のゲノムのリファレンス取得

リファレンス配列のfastaのみではなく、マッピングソフトのインデックスファイルや遺伝子情報ファイルも一緒に圧縮されて公開しています。
Copyright © Amelieff Corporation All Rights Reserved.

7

TRY!
RNA-seq解析：データ
データ取得 → クオリティコントロール → マッピング→発現定量• 酵母のゲノムのリファレンス配列を確認$ cd /home/ユーザ名/Desktop/amelieff/Scerevisiae$ llCopyright © Amelieff Corporation All Rights Reserved.

8

TRY!
RNA-seq解析：データ
データ取得 → クオリティコントロール → マッピング→発現定量• 酵母のゲノムのリファレンス配列を確認$ ll Bowtie2Index

今回はこちらのインデックスを使用します。
Copyright © Amelieff Corporation All Rights Reserved.

9RNA-seq解析：データ
データ取得 → クオリティコントロール → マッピング→発現定量• シーケンスデータ取得

http://trace.ddbj.nig.ac.jp/dra/index.html

DDBJのSequence Read Archive → Search

Copyright © Amelieff Corporation All Rights Reserved.

10RNA-seq解析：データ
データ取得 → クオリティコントロール → マッピング→発現定量• シーケンスデータ取得

Accessionに「 SRR518891 」と入力
Copyright © Amelieff Corporation All Rights Reserved.

→Search11

RNA-seq解析：データ
データ取得 → クオリティコントロール → マッピング→発現定量• シーケンスデータ取得

今回は1サンプルで実行しますが、発現比較する場合には複数サンプル必要で、replicateも多いほうが良いです。

ここからダウンロード

実験の詳細

NavigationエリアのExperiment→「SRX157933」をクリック

Copyright © Amelieff Corporation All Rights Reserved.

12RNA-seq解析：データ
データ取得 → クオリティコントロール → マッピング→発現定量• シーケンスデータ取得
他にも、シーケンサのプラットフォームやリード長などの情報も記載されています。

転写産物

Copyright © Amelieff Corporation All Rights Reserved.

13RNA-seq解析：データ
データ取得 → クオリティコントロール → マッピング→発現定量• シーケンスデータ取得（実行済み）
ダウンロードします。
$ wgetftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/fastq/SRA055/SRA055683/SRX157933/SRR518891_1.fastq.bz2$ wgetftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/fastq/SRA055/SRA055683/SRX157933/SRR518891_2.fastq.bz2

Copyright © Amelieff Corporation All Rights Reserved.

14

RNA-seq解析：データ
データ取得 → クオリティコントロール → マッピング→発現定量• シーケンスデータ取得（実行済み）
解凍して、先頭1000リードを抽出します。
$ bunzip2 SRR518891_1.fastq.bz2$ bunzip2 SRR518891_2.fastq.bz2$ head -4000 SRR518891_1.fastq > 1K_SRR518891_1.fastq$ head -4000 SRR518891_2.fastq > 1K_SRR518891_2.fastq

Copyright © Amelieff Corporation All Rights Reserved.

15

TRY!
RNA-seq解析：データ
データ取得 → クオリティコントロール → マッピング→発現定量• シーケンスデータを確認$ cd /home/ユーザ名/Desktop/amelieff$ llCopyright © Amelieff Corporation All Rights Reserved.

16

TRY!
RNA-seq解析：クオリティコントロール
データ取得 → クオリティコントロール → マッピング→発現定量• シーケンスデータのクオリティを確認
FastQCを実行します。
$ mkdir rnaseq$ fastqc -o rnaseq -f fastq 1K_SRR518891_1.fastq 1K_SRR518891_2.fastq

fastqc_report.htmlを、ウェブブラウザで開きます。
$ firefox rnaseq/1K_SRR518891_1_fastqc/fastqc_report.html$ firefox rnaseq/1K_SRR518891_2_fastqc/fastqc_report.html
以降の解析は、片側のリードのみ使用します。
Copyright © Amelieff Corporation All Rights Reserved.

17

TRY!
RNA-seq解析：クオリティコントロール
データ取得 → クオリティコントロール → マッピング→発現定量

PolyA/T tail が存在

最初の１塩基Copyright © Amelieff Corporation All Rights Reserved.

18

応用）

PolyA/T

tailの混入

RNA-seq（mRNA）では3'末端にPolyA/T tailがついている転写物をシーケンシングするため

Copyright © Amelieff Corporation All Rights Reserved.

19

TRY!
RNA-seq解析：クオリティコントロール
データ取得 → クオリティコントロール → マッピング→発現定量• 最初の１塩基を削除
fastx_trimmerの使用方法を確認する$ fastx_trimmer -h

Copyright © Amelieff Corporation All Rights Reserved.

20

TRY!
RNA-seq解析：クオリティコントロール
データ取得 → クオリティコントロール → マッピング→発現定量• 最初の１塩基を削除
2塩基目から使う$ cd rnaseq$ fastx_trimmer -f 2 -i ../1K_SRR518891_1.fastq -o 1K_SRR518891_1_s.fastq -Q33

Copyright © Amelieff Corporation All Rights Reserved.

21

TRY!
RNA-seq解析：クオリティコントロール
データ取得 → クオリティコントロール → マッピング→発現定量• PolyA/T tailを除去
3'端にAを5連続以上含むリード数がどのくらいあるか調べる$ grep "AAAAA$" 1K_SRR518891_1_s.fastq | wc -l
39
fastx_clipperの使用方法を確認する$ fastx_clipper –hCopyright © Amelieff Corporation All Rights Reserved.

22

TRY!
RNA-seq解析：クオリティコントロール
データ取得 → クオリティコントロール → マッピング→発現定量• PolyA/T tailを除去
PolyA/T tailを除去する$ fastx_clipper -a AAAAA -i 1K_SRR518891_1_s.fastq-o 1K_SRR518891_1_s_notail.fastq -Q 33

Prinseqなど、各リードのPolyA/Tの数に合わせて除去するソフトもあります。

Copyright © Amelieff Corporation All Rights Reserved.

23

TRY!
RNA-seq解析：クオリティコントロール
データ取得 → クオリティコントロール → マッピング→発現定量• クオリティの低いリードを除外
3’末端からクオリティ20未満の塩基をトリミングし、長さが30塩基未満になったリードを破棄する。
さらに、80%以上の塩基がクオリティー20以上のリードのみを抽出する。
$ fastq_quality_trimmer -t 20 -l 30 -Q 33-i 1K_SRR518891_1_s_notail.fastq | fastq_quality_filter-q 20 -p 80 -Q 33 -o 1K_SRR518891_1_clean.fastq

Copyright © Amelieff Corporation All Rights Reserved.

24

TRY!
RNA-seq解析：クオリティコントロール
データ取得 → クオリティコントロール → マッピング→発現定量• クリーニング結果の確認
FastQCを実行します。
$ fastqc -f fastq 1K_SRR518891_1_clean.fastq

fastqc_report.htmlを、ウェブブラウザで開きます。
$ firefox 1K_SRR518891_1_clean_fastqc/fastqc_report.html

クリーニング前後のリード数と、クオリティの変化を確認してください。
Copyright © Amelieff Corporation All Rights Reserved.

25

TRY!
RNA-seq解析：クオリティコントロール
データ取得 → クオリティコントロール → マッピング→発現定量

クリーニング前

クリーニング後

サンプルや調整方法、シーケンサの特徴にあわせてクリーニング項目や条件を工夫しています。
Copyright © Amelieff Corporation All Rights Reserved.

26

応用）

アダプタ配列を除外するソフト① cutadapt② FastX-Toolkit（fastxclipper）③ tagcleaner

指定した配列はどのソフトでも除ける。
fastx_clipperは、部分配列もかなり除けたが、リード数も1/10以下に減るためアダプタ以外の配列も除いている可能性がある。
tagcleanerは、一度に1アダプタ配列しか指定できない。
$ cutadapt-b TCTCGTATGCCGTCTTC -b CTACAGTCCGACGA-m 10 -n 2 $FILE.fastq 1> $FILE_cutadapt.fastq※オプションmnOe

これより短くなったものは破棄
同リードへのアダプタ出現回数
マッチ領域の最少長「エラー塩基数/マッチ領域長」の最大
Copyright © Amelieff Corporation All Rights Reserved.

27

TRY!
RNA-seq解析：マッピング
データ取得 → クオリティコントロール → マッピング→発現定量• TopHatの使い方を確認$ tophat

スプライシングを考慮して、マッピングするため、既知の遺伝子情報を使用することもできます。

Copyright © Amelieff Corporation All Rights Reserved.

28

TRY!
RNA-seq解析：マッピング
データ取得 → クオリティコントロール → マッピング→発現定量• マッピング※今回はリード数が少ないため、マッピング基準を緩めています。

$ tophat -o 1K_SRR518891 -g 3 -N 10 --read-edit-dist 10--read-gap-length 10/home/ユーザ名/Desktop/amelieff/Scerevisiae/Bowtie2Index/genome
1K_SRR518891_1_clean.fastq

BAMとインデックス、BEDなどが作成されます。

$ ls 1K_SRR518891-N/--read-mismatches--read-edit-dist--read-gap-length

Final read alignments having more than these many mismatches are discarded.
Final read alignments having more than these many edit distance are discarded.
Final read alignments having more than these many total length of gaps are discarded.
Copyright © Amelieff Corporation All Rights Reserved.

29

TRY!
RNA-seq解析：マッピング
データ取得 → クオリティコントロール → マッピング→発現定量• マッピング率$ less 1K_SRR518891/align_summary.txt

マッピング率Copyright © Amelieff Corporation All Rights Reserved.

30

ポ イ ン ト ） To p H a t の ア ル ゴ リ ズ ム
1. リードをペアエンドでリファレンスにマッピングする。
2. マッピングできなかったリードを断片化して、リファレンスにマッピングする。
3. マッピング結果をもとに、転写構造をアセンブリングする。

http://en.wikipedia.org/wiki/File:RNA-seq-alignment.png

http://www.ncbi.nlm.nih.gov/pubmed/19289445

Copyright © Amelieff Corporation. All Rights Reserved.

31

TRY!
RNA-seq解析：マッピング
データ取得 → クオリティコントロール → マッピング→発現定量• マッピング結果の可視化$ samtools index accepted_hits.bam$ igv.shaccepted_hits.bamをIGVで表示してください。

Copyright © Amelieff Corporation All Rights Reserved.

32

TRY!
RNA-seq解析：マッピング
データ取得 → クオリティコントロール → マッピング→発現定量• マッピング結果の可視化

Positionの例：
III:177,350-177,500

Copyright © Amelieff Corporation All Rights Reserved.

33

ポイント）
遺伝子の発現量 ≠ 遺伝子上にマップされたリード数

長い遺伝子ほど、マップされるリードは多くなる（遺伝子間のバイアス）
サンプル量の多いランほど、マップされるリードは多くなる（ラン間のバイアス）
これらのバイアスを補正してから発現量を比較する必要があります

・発現量としてよく使われる指標
RPKM（Reads Per Kilobase per Million mapped reads）
FPKM（Fragments Per Kilobase of exon per Million mapped fragments）
どちらも、発現量をエクソン長と全マッピング数で補正した値

1,000,000
1,000
FPKM = raw counts ××
all reads gene length
Copyright © Amelieff Corporation All Rights Reserved.

34

TRY!
RNA-seq解析：発現定量
データ取得 → クオリティコントロール → マッピング→発現定量• Cufflinksの使い方を確認$ cufflinks

アセンブルのガイドとして既知の遺伝子情報を使用することもできます。
Copyright © Amelieff Corporation All Rights Reserved.

35

TRY!
RNA-seq解析：発現定量
データ取得 → クオリティコントロール → マッピング→発現定量• 発現量を計算※今回はリード数が少ないため、検出基準を緩めています。

$ cufflinks --min-frags-per-transfrag 2-o 1K_SRR518891 1K_SRR518891/accepted_hits.bam

fpkm_trackingファイルが作成されます。

$ ll –h 1K_SRR518891--min-frags-per-transfrag minimum number of fragments needed for new transfrags

Copyright © Amelieff Corporation All Rights Reserved.

36

TRY!
RNA-seq解析：発現定量
データ取得 → クオリティコントロール → マッピング→発現定量• 発現量を計算$ less 1K_SRR518891/genes.fpkm_tracking

4列目がGene ID、10列目がFPKMです。

Copyright © Amelieff Corporation All Rights Reserved.

37応用）サンプル間比較
サンプルごとに発現量を計算したあと、サンプルごとに発現している遺伝子が違うため、比較の基準とする遺伝子リストを作成します。
$ cuffmerge -o COMPARE -g genes.gtf -s genome.fa
transcript.gtf.txt
Group1/S1/transcript.gtf
Group1/S2/transcript.gtf
Group2/S3/transcript.gtf
Group2/S4/transcript.gtf

transcript.gtf.txt
各サンプルのcufflinks結果を羅列したファイル

発現量を比較します$ cuffdiff -o COMPARE -L Group1,Group2 genes.gtf
Group1/S1/accepted_hits.bam,Group1/S2/accepted_hits.bam
Group2/S3/accepted_hits.bam,Group2/S4/accepted_hits.bam
Copyright © Amelieff Corporation All Rights Reserved.

38本講義の内容• Reseq解析• RNA-seq解析

公開データ取得

公開データ取得↓↓

クオリティコントロール

クオリティコントロール

↓

↓

マッピング

マッピング

↓

↓

変異検出

発現定量

Copyright © Amelieff Corporation All Rights Reserved.

39

ご清聴ありがとうございました。

Copyright © Amelieff Corporation All Rights Reserved.

40