機能ゲノム学（第6回）

Size: px

Start display at page:

Download "機能ゲノム学（第6回）"

まとものあき
4 years ago
Views:

1 R でトランスクリプトーム解析東京大学大学院農学生命科学研究科アグリバイオインフォマティクス教育研究ユニット門田幸二 ( かどたこうじ ) kadota@iu.a.u-tokyo.ac.jp 1

生命情報科学研究センター (CBRC) 産総研特別研究員 2003/11/1~ 放医研先端遺伝子発現研究センター研究員

2 自己紹介 2002 年 3 月東京大学大学院農学生命科学研究科博士課程修了学位論文 : cdna マイクロアレイを用いた遺伝子発現解析手法の開発 ( 指導教官 : 清水謙多郎教授 ) 2002/4/1~ 産総研生命情報科学研究センター (CBRC) 産総研特別研究員 2003/11/1~ 放医研先端遺伝子発現研究センター研究員 2005/2/16~ 東京大学大学院農学生命科学研究科特任助手 2007/4/1~ 現在東京大学大学院農学生命科学研究科特任助教アグリバイオインフォマティクスプログラム 2

3 様々な Motivation ~ の原因遺伝子 ( ガン関連遺伝子とか ) を同定したい FASTQ 以降の一通りの解析ができるようになりたい (Windows の )R でできることとできないことモデル生物と非モデル生物の解析戦略の違い倍率変化で解析 vs. 分布を使って解析いろんな R パッケージがあるけれど RNA-seq で二つのサンプルを比較し発現変動遺伝子同定までを行うまでの流れを一通り紹介 A 群腎臓正常組織 wildtype B 群肝臓腫瘍組織 mutant 3

4 データ解析のスタート地点 NGS から得られた FASTQ 形式ファイルデータ取得完了! A1.fq A2.fq B1.fq B2.fq なんじゃこの変な記号は! 何をどうすれば... 4

5 比較トランスクリプトーム解析の流れ複数の FASTQ ファイル Linux マシンリファレンス配列の作成 A1.fq, A2.fq, B1.fq, B2.fq 複数サンプルの混合アセンブルにより RefSeq のような転写物配列集合 (multi-fasta ファイル ) を得るイメージマッピングどの転写物にどれだけの数のリードがマップされたかといういわゆる遺伝子発現行列を得るイメージデータ解析発現変動遺伝子のリストアップや作図などすべてが (Windows の )R でできるわけではありません 5

6 Linux マシン使用部分の解決策自前で大容量メモリ計算サーバ (Linux) を購入し必要なソフトのインストールからスタート特徴 : 難易度は高いが思い通りの解析が可能 Linux サーバをもつバイオインフォ系の人にお願いする特徴 : 気軽に頼める知り合いがいればいいがその人次第 DDBJ Read Annotation Pipeline を利用特徴 : 一番お手軽な選択肢だがサポートしているプログラムのみデータ登録が前提?! だが手取り足取り丁寧に教えてくれるので個人的にはこちらを推奨自分の負担も減るし 6

7 7

8 比較トランスクリプトーム解析の流れ複数の FASTQ ファイル Linux マシンリファレンス配列の作成クオリティチェックアセンブル結果 (multi-fasta) ファイルから平均長やトータルの長さなどの基本情報を抽出マッピングマッピング結果 (BED 形式 ) ファイルを入力として転写物ごとのマップされたリード数をカウントデータ解析発現変動遺伝子のリストアップや作図など大規模計算部分以外は一通りできます 8

9 参考ウェブページ 9

10 比較トランスクリプトーム解析の流れ複数の FASTQ ファイルクオリティチェックリファレンス配列の作成アセンブル結果 (multi-fasta) ファイルから平均長やトータルの長さなどの基本情報を抽出マッピングマッピング結果 (BED 形式 ) ファイルを入力として転写物ごとのマップされたリード数をカウントデータ解析発現変動遺伝子のリストアップや作図など 10

11 データのクオリティチェックデスクトップにあるフォルダ中に SRR fastq というファイルが存在するという前提 11

12 R の起動デスクトップにあるフォルダ中のファイルを解析 12

13 作業ディレクトリ (= フォルダ ) の変更

14 getwd() と打ち込んで確認 14

15 コピー & ペースト ( こぴぺ ) 一連のコマンド群をコピーして 2R Console 画面上でペースト 15

16 html レポートが作成される R ではうまくいくが R ではエラーがでます 16

17 html レポートが作成される (R では ) ポジションごとのクオリティスコアなどの情報が得られます 17

: 配列情報三行目 : + からはじまる一行 ( の description 行 ) 四行目 : クオリティ情報 @SEQ_ID

18 FASTQ 形式 ( と FASTA 形式 ) FASTA 形式 > ではじまる一行の description 行と配列情報からなる形式 NGS の read 長は短いので実質的に一つのリードを二行で表現 >SEQ_ID GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT FASTQ 形式一行目ではじまる一行の description 行二行目 : 配列情報三行目 : + からはじまる一行 ( の description 行 ) 四行目 : GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT +!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC

19 塩基配列のクオリティ情報といえば Phred スコア Phred というベースコールプログラムから得られる Quality Value(QV 値 ) のことなぜ FASTQ 形式では Phred スコアそのものでクオリティ情報を表現しないの? 19

20 理由 :( 容量 ) 節約のため Cock et al., Nucleic Acids Res., 38: , 2010 FASTQ 形式中のクオリティ情報部分 Phred スコア (QUAL 形式 ) Phred スコアが X の場合 ASCII (X+33) に対応する文字コードを割り当てる 20

21 R の現実 (R で ) 塩基配列解析のウェブページは常に最新の情報が記載されているわけではない昔のバージョンだとうまくいくが最新版だとエラーがでることもある R のバージョンを上げてから昔作ったスクリプトを実行しようとするとそんな関数ないと文句を言われたよく調べてみると関数名が変わっていた例 :DESeq パッケージ中の estimatevariancefunctions estimatedispersions DEGseq のスクリプトがうまく動かなくなっている R ( ) R ( ) での体験談 21

22 R の現実と対処法ぐぐる ( qrqc error などで検索) Rのバージョンを最新版 ( または古いもの ) に変更 sessioninfo() で現在利用している R のバージョンを確認手元にある印刷物のマニュアルは古いものでは?! 22

23 R の昔のバージョンのインストール (R で ) マイクロアレイデータ解析のウェブページのほうです 23

24 比較トランスクリプトーム解析の流れ複数の FASTQ ファイルクオリティチェックリファレンス配列の作成アセンブル結果 (multi-fasta) ファイルから平均長やトータルの長さなどの基本情報を抽出マッピングマッピング結果 (BED 形式 ) ファイルを入力として転写物ごとのマップされたリード数をカウントデータ解析発現変動遺伝子のリストアップや作図など 24

25 リファレンス配列について Q: マッピングに使うリファレンス配列は? A: 好きなものを使ってくださいゲノム配列でもトランスクリプトーム配列でも結構です Q: どこから取得できるんですか? A: UCSC Sequence and Annotation Downloads などから取得できます ( アノテーション情報も ) 以下はほんの一例ヒト全ゲノム配列の場合 ftp://genome-ftp.cse.ucsc.edu/goldenpath/hg19/bigzips/hg19.2bit ヒトトランスクリプトーム配列 (RefSeq mrna) の場合 ftp://genome-ftp.cse.ucsc.edu/goldenpath/hg19/bigzips/refmrna.fa.gz ヒトアノテーション情報の場合

gov/refseq/h_sapiens/ mrna_prot/human.rna.fna.

26 こんな感じヒトトランスクリプトーム配列ヒトゲノム配列 1-22 番染色体 +X+Y 約 6200 万行のファイル約 3GB のサイズ 46,XXX mrna sequences ftp://ftp.ncbi.nih.gov/refseq/h_sapiens/ mrna_prot/human.rna.fna.gz アノテーション情報ファイル (refflat 形式 ) symbol 染色体名転写開始位置 name CDS start Exon 数 strand 転写終結位置 CDS end 26

27 非モデル生物の場合手持ちの RNA-Seq データのみ 2010 年頃から提供されはじめた de novo transcriptome assembly 用のプログラム (Trinity や Trans-ABySS など ; もちろん Linux 用です ) を利用すればトランスクリプトームの配列セット (RefSeq のようなイメージ ) を得ることができます入力 :RNA-Seq データ出力 : コンティグ ( 転写物配列 ) >read1 GGGGTTCAAAGCAGTATCGATCAAATAGTA >read2 GTTCAAAGCAGTATCGATCAAATAGTAAAT >read3 ACGATGCAGCCTTAACGATGGTCCACAATT >read4 >contig1(transcript1) GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAA CTCACAGTTTGGAGCTTATCAGTCAA >contig2(transcript2) ACGATGCAGCCTTAACGATGGTCCACAATTATCGGGAATCA >contig3(transcript3) 27

28 非モデル生物の場合手持ちの RNA-Seq データのみでアセンブルを行う場合は paired-end のデータが基本ペアードエンド法断片配列の両末端が数百塩基以内の対の二種類の配列が得られる A1: A1_1.fq A1_2.fq A2: A2_1.fq A2_2.fq シングルエンド法約塩基 B1: B1_1.fq B1_2.fq B2: B2_1.fq B2_2.fq read_1.fq read_2.fq 二つのファイルを入力として比較対象サンプルを混合したデータのアセンブルリファレンストランスクリプトーム配列の取得 28

29 非モデル生物の場合 Trinity (Grabherr et al., Nat. Biotechnol., 2011) 実行プログラムの一例 nohup Trinity.pl --seqtype fq --left read_1.fq --right read_2.fq -- run_butterfly --bflyheapspace 180G --CPU 2 --bfly_opts "-V stderr" --run_allpathslg_error_correction --output hoge & アセンブルが終了すると hoge というディレクトリ中に Trinity.fasta という multi-fasta ファイルが得られる 29

30 比較トランスクリプトーム解析の流れ複数の FASTQ ファイルクオリティチェックリファレンス配列の作成アセンブル結果 (multi-fasta) ファイルから平均長やトータルの長さなどの基本情報を抽出マッピングマッピング結果 (BED 形式 ) ファイルを入力として転写物ごとのマップされたリード数をカウントデータ解析発現変動遺伝子のリストアップや作図など 30

31 multi-fasta 形式ファイルからの情報抽出一連のコマンド群をコピーして 2R Console 画面上でペースト 31

32 multi-fasta 形式ファイルからの情報抽出 1 出力ファイル名として指定したもの (hoge.txt) がフォルダ中に作成される ( はず ) 32

33 練習中にある practice1.txt 中の記述を変更して Trinity.fasta ファイルに対して同様の解析を行い結果を hoge1.txt に出力せよ 33

34 multi-fasta 形式ファイルからの情報抽出 2 配列ごとの GC 含量を計算したいとき 34

35 練習中にある practice2.txt 中の記述を変更して Trinity.fasta ファイルに対して同様の解析を行い結果を hoge2.txt に出力せよ 35

36 multi-fasta 形式ファイルからの情報抽出 3 Trinity.fasta に対して 600bp 以上のもののみ抽出してみよう 36

37 multi-fasta 形式ファイルからの情報抽出 4 FPKM 値 : 配列長補正ずみの発現量に相当する値 Trinity の最近のバージョンでは FPKM 値を出力しなくなったようですね 37

38 multi-fasta 形式ファイルからの情報抽出 4 Trinity.fasta FPKM 値をもとにサンプル内の転写物間の発現レベルの大小を議論可能サンプル間の比較には使えない ( といわれている ) 38

39 比較トランスクリプトーム解析の流れ複数の FASTQ ファイルクオリティチェックリファレンス配列の作成アセンブル結果 (multi-fasta) ファイルから平均長やトータルの長さなどの基本情報を抽出マッピングマッピング結果 (BED 形式 ) ファイルを入力として転写物ごとのマップされたリード数をカウントデータ解析発現変動遺伝子のリストアップや作図など 39

40 マッピングの基本的なイメージ基本的なマッピングプログラム (bowtie など ) を用いた場合リファレンス配列 : ゲノム count T1 サンプルの RNA-Seq データ mapping 遺伝子 1 遺伝子 2 遺伝子 3 遺伝子 4 リファレンス配列 : トランスクリプトーム count 遺伝子 1 遺伝子 2 遺伝子 3 遺伝子 4 40

41 解析の目的に応じてプログラムを使い分けゲノム配列既知で遺伝子構造を知りたいような場合 Cufflinks (Trapnell et al., Nat. Biotechnol., 2010) ARTADE2 (Kawaguchi et al., Bioinformatics, 2012) ( ゲノム配列の有無に関わらず )RefSeq のようなトランスクリプトーム配列にマッピングして比較トランスクリプトーム解析を行いたい場合 Bowtie (Langmead et al., Genome Biol., 2009) BWA (Li and Durbin, Bioinformatics, 2009) Cufflinks 周辺は専門外です 41

42 マッピング結果の出力ファイル形式 ( ゲノム配列の場合 ) どの染色体上のどの位置に ( どのリードが ) マッピングされたかあるいは ( トランスクリプトーム配列の場合 ) どの転写物配列上のどの位置に ( どのリードが ) マッピングされたかを表すファイル形式 ( フォーマット ) は複数あります : BED (Browser Extensible Data) format BEDtools (Quinlan et al., Bioinformatics, 26: , 2010) GFF (General Feature Format) format SAM (Sequence Alignment/Map) format SAMtools (Li et al., Bioinformatics, 25: , 2009) マッピング結果ファイルからどうやって転写物ごとのマップされたリード数をカウントするのか? 42

43 BED 形式 43

44 BED 形式あるトランスクリプトーム配列 (RefSeq) にマップした結果転写物 ID Start End 転写物 ID ごとの出現数 = マップされたリード数 44

45 比較トランスクリプトーム解析の流れ複数の FASTQ ファイルクオリティチェックリファレンス配列の作成アセンブル結果 (multi-fasta) ファイルから平均長やトータルの長さなどの基本情報を抽出マッピングマッピング結果 (BED 形式 ) ファイルを入力として転写物ごとのマップされたリード数をカウントデータ解析発現変動遺伝子のリストアップや作図など 45

46 R を用いてコピペでマップされたリード数情報を得ることができます output3.txt 46

47 データ解析の前に研究目的 ( と手持ちのデータ ) をおさらい一つのサンプル内でどの転写物 (or 遺伝子 ) の発現レベルが高いか低いかを調べたい場合 RPKM や FPKM などの転写物の長さを考慮して正規化されたデータで解析トータルのリード数を補正する必要はないがやってもよい遺伝子間の発現レベルの大小関係を調べたいだけなので解析データを定数倍したところで何ら影響を与えないからサンプル間比較 (sample A vs. B など ) で発現変動遺伝子 ( Differentially Expressed Genes; DEGs) を調べたい場合トータルのリード数を補正したデータで解析正確にはサンプル間で発現変動していない遺伝子 (non-degs) ができるだけ発現変動していないと判定されるように正規化したデータ既存の R パッケージを用いて解析を行う場合には ( 整数値のみからなる ) 生のリードカウントデータを入力とし内部的に上記正規化を行う研究目的によってやっていい正規化とやってはいけない ( と言われている ) 正規化がある 47

48 比較トランスクリプトーム解析の流れ複数の FASTQ ファイルクオリティチェックリファレンス配列の作成アセンブル結果 (multi-fasta) ファイルから平均長やトータルの長さなどの基本情報を抽出マッピングマッピング結果 (BED 形式 ) ファイルを入力として転写物ごとのマップされたリード数をカウントデータ解析発現変動遺伝子のリストアップや作図など 48

49 二群間比較用 R パッケージ DEGseq (Wang et al., Bioinformatics, 26: , 2010) ポワソン分布 (variance = mean) を仮定しているためばらつきを過少評価 edger (Robinson et al., Bioinformatics, 26: , 2010) 正規化法 :TMM 法負の二項分布 (variance > mean) を仮定 mean のみのパラメータを用いて現実のばらつきを表現 DESeq (Anders and Huber, Genome Biol., 11: R106, 2010) 正規化法 :RLE 法 (relative log expression) edger のモデルをさらに拡張 ( しているらしい ) bayseq (Hardcastle and Kelly, BMC Bioinformatics, 11:422, 2010) 正規化法 :RPM ( たぶん ) 配列の長さ情報を与えるオプションがあるデータセット中に占める DEG の割合 (P DEG ) を一意に返す NBPSeq (Di et al., SAGMB, 10:24, 2011) 入力 : 生のリードカウントからなる遺伝子発現行列出力 : 遺伝子ごとの発現変動の度合い (p 値など ) 49

50 理想的な実験デザイン ( 二群間比較 ) サンプル A vs. B の比較 (Kidney vs. Liver; 腎臓 vs. 肝臓 ) 生のリードカウントのデータ ( 整数値 ) Biological replicates のデータ生物学的なばらつき ( 個体間の違い ) を考慮すべし A1: ある生物の腎臓 A2: 同じ生物種の別個体の腎臓 A3: 同じ生物種のさらに別個体の腎臓 B1: ある生物の肝臓 B2: 同じ生物種の別個体の肝臓 50

51 分布の話例題 :Marioni et al., Genome Res., 18: , 2008 のデータ ( の一部 ) kidney( 腎臓 ) liver( 肝臓 ) Technical replicates のデータサンプル内の技術的なばらつき ( 例 : レーン間の違い ) の度合いを調べるためのデータでありこのようなデータで二群間比較し発現変動遺伝子がどの程度あるかといった数に関する議論は無意味解析例 : アリエナイ?! 数 (50% とか ) が発現変動遺伝子として検出される理由 :Biological variation > Technical variation 51

52 分布の話例題 :Marioni et al., Genome Res., 18: , 2008 のデータ ( の一部 ) kidney( 腎臓 ) RPM 正規化 1,000,000 12,685 1,804,

分布の話例題 :Marioni et al., Genome Res., 18: 1509-1517, 2008 のデータ ( の一部 ) kidney( 腎臓 ) adjusted R-squared: 0.

53 分布の話例題 :Marioni et al., Genome Res., 18: , 2008 のデータ ( の一部 ) kidney( 腎臓 ) adjusted R-squared: y = x y = a + bx Technical replicates のデータは : ( 遺伝子の )VARIANCE はその MEAN で説明可能である VARIANCE MEAN ポアソン分布に従うポアソンモデルが適用可能 53

54 分布の話生物アイコン ( 例題 :Cumbie et al., PLoS ONE, 6: e25279, 2011 のデータ ( の一部 ) Arabidopsis( シロイヌナズナ ) adjusted R-squared: y = a + bx y = x Biological replicates のデータは : VARIANCE > MEAN 負の二項 (NB) 分布に従う NB モデルが適用可能 54

55 なぜ沢山の方法が存在しうるのか? DEGseq (Wang et al., Bioinformatics, 26: , 2010) ポワソン分布 (variance = mean) を仮定しているためばらつきを過少評価 edger (Robinson et al., Bioinformatics, 26: , 2010) 正規化法 :TMM 法負の二項分布 (variance > mean) を仮定 DESeq (Anders and Huber, Genome Biol., 11: R106, 2010) 正規化法 :RLE 法 (relative log expression) edgerのモデルをさらに拡張 ( しているらしい ) bayseq (Hardcastle and Kelly, BMC Bioinformatics, 11:422, 2010) 正規化法 :RPM ( たぶん ) 配列の長さ情報を与えるオプションがあるデータセット中に占めるDEGの割合 (P DEG ) を一意に返す NBPSeq (Di et al., SAGMB, 10:24, 2011) VAR (1 ) VAR (1 ) Ans. Variance と Mean の関係を表現する手段が沢山あるから VAR (1 1 ) VAR 55

56 edger を使ってみる例題 :Marioni et al., Genome Res., 18: , 2008 のデータ kidney( 腎臓 ) liver( 肝臓 ) ファイル名 :SupplementaryTable2_changed.txt 内容 :A 群が最初の 5 列 B 群が残りの 5 列のデータ解析結果を hoge2.txt という名前でファイルに出力したい 56

57 edger を使ってみるファイル名 :SupplementaryTable2_changed.txt 内容 :A 群が最初の 5 列 B 群が残りの 5 列のデータ解析結果を hoge2.txt という名前でファイルに出力したい 57

58 edger を使ってみる R 上でスクリプトをコピペ! ( エラーメッセージが出ていなければ hoge2.txt というファイルができているはず ) 58

59 edger を使ってみる一番右側の数値が False Discovery Rate (FDR) この列 (O 列 ) で昇順にソートすれば任意の閾値を満たす遺伝子数がわかる 19,785 個が FDR < 0.01 を満たす 21,291 個が FDR < 0.05 を満たす 59

60 edger を使ってみる Top-ranked gene の生リードカウントを眺めても確かに発現変動 (Kidney << Liver) していることが分かる 60

61 edger を使ってみる M-A plot を描画 (FDR < 0.01 を満たすものを赤色で表示 ) hoge2.png 19,794 個 ( 全遺伝子数のうち約 62% が FDR < 0.01 を満たす ) 61

62 edger を使ってみる M-A plot を描画 (2 倍以上発現変動しているものを赤色で表示 ) hoge2.png 個 ( 全遺伝子数のうち約 37% が 2 倍以上発現変動している ) このやり方はダメなんです 62

63 倍率変化がだめな理由をデモ例題 :Marioni et al., Genome Res., 18: , 2008 のデータ kidney( 腎臓 ) liver( 肝臓 ) 発現変動遺伝子がないデータで二群間比較をしてみる A 群 B 群 63

64 倍率変化がだめな理由をデモ例題 :Marioni et al., Genome Res., 18: , 2008 のデータ ( の一部 ) (A1, A2) vs. (A3, A4) の二群間比較結果 edger で FDR < 0.01 を満たすものは 0 個 (edger で )2 倍以上発現変動しているものは 3814 個 Rcode_edgeR_tech_rep_fdr001.txt Rcode_edgeR_tech_rep_fc2.txt 低発現領域で log 比が大きくなる現象をうまくモデル化することが重要 64

65 Top 400 Top 2000 低い全体的な発現レベル高いこんな感じでランキングすることが重要です 65

66 26,221 genes Biological replicates の 3 vs. 3 サンプル例題 :Cumbie et al., PLoS ONE, 6: e25279, 2011 の Arabidopsis データ A 群 B 群 data_arab.txt オリジナルは AT4G32850 のものが重複して存在していたため行目のデータを予め除去している 66

67 edger を default の手順 (edger/default) で実行 67

68 edger を default の手順 (edger/default) で実行 B 群で高発現 A 群で高発現 68

69 サンプル間クラスタリングも重要です 69

70 サンプル間クラスタリングも重要ですデータ中に発現変動遺伝子がありそうかどうかはクラスタリング結果を眺めるだけでかなりわかる 70

71 東大生以外の方も受講可能です ( 来年度もやります ) 71

72 謝辞共同研究者清水謙多郎先生 ( 東京大学 ) 嶋田透先生 ( 東京大学 ) 西山智明先生 ( 金沢大学 ) 勝間進先生 ( 東京大学 ) 河岡慎平博士 ( 東京大学 ) 末次克行先生 ( 農業生物資源研究所 ) 上樂明也先生 ( 農業生物資源研究所 ) グラント若手研究 (B)(H21-23 年度 ): マイクロアレイ解析の再現性感度特異度を飛躍的に向上させるデータ解析手法の開発 ( 代表 ) 新学術領域研究 ( 研究領域提案型 )(H22 年度 -): 非モデル生物におけるゲノム解析法の確立 ( 分担 ; 研究代表者 : 西山智明 ) 72

機能ゲノム学（第6回）

機能ゲノム学（第6回） RNA-Seq データ解析リテラシー東京大学大学院農学生命科学研究科アグリバイオインフォマティクス教育研究ユニット門田幸二 ( かどたこうじ ) http://www.iu.a.u-tokyo.ac.jp/~kadota/ kadota@iu.a.u-tokyo.ac.jp 1 2009 年ごろの私次世代シーケンサー (NGS) 解析についての認識単に短い塩基配列が沢山あるだけでしょ得られる配列データって