ゲノム情報解析基礎～ Rで塩基配列解析～

Size: px

Start display at page:

Download "ゲノム情報解析基礎～ Rで塩基配列解析～"

しらんおとべ
5 years ago
Views:

1 ネット接続できないヒトもダブルクリックでローカルに r_seq.html を起動可能です実習はデスクトップ上にある hoge フォルダの中身が以下の状態を想定して行います (R で ) 塩基配列解析の利用法 : GC 含量計算から発現変動解析まで東京大学大学院農学生命科学研究科アグリバイオインフォマティクス教育研究プログラム門田幸二 ( かどたこうじ ) kadota@iu.a.u-tokyo.ac.jp 1

2 参考ウェブページ (R で ) 塩基配列解析の利用法の紹介 2

3 前提条件以下の手順通りに R および必要なパッケージのインストールが完了しているという前提です 3

産総研特別研究員マイクロアレイ解析手法開発 2003/11/1~ 放医研先端遺伝子発現研究センター研究員一次元電気泳動波形解析手法開発

4 自己紹介 2002 年 3 月東京大学大学院農学生命科学研究科応用生命工学専攻博士課程修了学位論文 : cdna マイクロアレイを用いた遺伝子発現解析手法の開発 ( 指導教官 : 清水謙多郎教授 ) 2002/4/1~ 産総研生命情報科学研究センター (CBRC) 産総研特別研究員マイクロアレイ解析手法開発 2003/11/1~ 放医研先端遺伝子発現研究センター研究員一次元電気泳動波形解析手法開発 2005/2/16~ 東京大学大学院農学生命科学研究科アグリバイオインフォマティクスプログラムマイクロアレイ解析手法開発 RNA-seq データ解析手法開発 ( トランスクリプトーム解析周辺の ) 手法開発系のヒトです 4

5 講義風景 ( 平成 26 年度 ) 5

6 Contents R でゲノム解析シロイヌナズナゲノムの GC 含量計算 multi-fastaファイルの読み込み関数やオプションの利用法パッケージの説明 R でトランスクリプトーム解析シロイヌナズナの RNA-seq データを一通り解析公共 DBからの生データ取得マッピングおよびカウントデータ取得サンプル間クラスタリング発現変動遺伝子 (DEG) 検出 6

7 植物グローバルなので例 : シロイヌナズナ (Arabidopsis thaliana) ゲノム配列決定 (chr1-5, chrc, chrm) 1 番染色体 :Theologis et al., Nature, 408: , 番染色体 :Lin et al., Nature, 402: , 番染色体 :Salanoubat et al., Nature, 408: , 2000 トランスクリプトーム配列 (cdna 配列 ) 決定アノテーション :Seki et al., Science, 296: , 2002 まとめサイト The Arabidopsis Information Resource (TAIR) Lamesch et al., Nucleic Acids Res., 40: D , R でゲノム解析が可能です 7

8 アノテーションファイル?! TAIR10 のゲノム配列ファイル (TAIR10_chr_all.fas) はここからダウンロードしました 8

9 アノテーションファイル?! 参考 TAIR ウェブインターフェースから取得する際のイメージ 9

10 参考 TAIR10 のゲノム配列ファイル (TAIR10_chr_all.fas) はこんな感じです 10

11 multi-fasta ファイル?! > のヘッダ行塩基またはアミノ酸配列が複数 (multi) 個からなるファイルのこと 11

12 R で配列長と GC 含量計算最新のゲノム配列を読み込んで GC 含量を計算することができます TAIR10_chr_all.fas は hoge フォルダ中にあり 12

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

14 作業ディレクトリの変更 Windows(C:) となっている場合もあるが気にしないこの部分はひとそれぞれ

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

基本はコピペ 1 2013 年 7 月以降のリニューアルでコードのコピーがやりずらくなっています CTRL と ALT

16 基本はコピペ年 7 月以降のリニューアルでコードのコピーがやりずらくなっています CTRL と ALT キーを押しながらコードの枠内で左クリックすると全選択できます 2 1 一連のコマンド群をコピーして 2R Console 画面上でペースト 16

17 実行結果実行前の hoge フォルダ実行後の hoge フォルダ出力 : hoge4.txt 17

18 ありがちなミス 1 作業ディレクトリの変更を忘れているため入力ファイルの読み込み段階でエラーとなる 18

19 ありがちなミス 2 必要な入力ファイルが作業ディレクトリ中に存在しない 19

20 ありがちなミス 3 出力予定のファイル名と同じものを別のプログラムで開いているため最後の write.table 関数のところでエラーが出る 20

21 ありがちなミス 4 実行スクリプトをコピーする際最後の行のところで改行を含ませずに R Console 画面上でペーストしたため最後のコマンドが実行されない ( 出力ファイルが生成されない ) 21

22 R で配列長と GC 含量計算出力 : hoge4.txt 原著論文中の数値ちゃんと似た結果が得られています 22

23 詳細を説明出力 : hoge4.txt 入力と出力ファイル名を指定しているところ 23

24 詳細を説明 GC 含量計算をしたいときには Biostrings というパッケージを読み込む必要がありますこの作業を行っておかないと例えば multi-fasta ファイルを読み込むための readdnastringset 関数を利用できません 24

25 詳細を説明 readdnastringset 関数を用いて 1in_f で指定した入力ファイルを 2fasta 形式で読み込んだ結果を 3fasta というオブジェクト名で格納しています

26 詳細を説明 fasta オブジェクトは確かに multi-fasta ファイル中の情報を適切に読み込めている 26

27 色についての説明 27

28 色についての説明単に確認しているだけなので灰色になっている 28

29 浮かんでくる疑問 FASTA 形式以外にどんな形式を読み込めるのだろう?FASTQ 形式は読み込めるの? DNA 配列以外にどんな配列を読み込めるのだろう? アミノ酸配列は? 29

30 関数の使用法について? 関数名で使用法を記したウェブページが開くページの下のほうに ( 大抵の場合 ) 使用例が掲載されている使用法既知の関数のマニュアルをいくつか読んで慣れておく 30

31 関数の使用法について FASTQ ファイルは読み込めそうだ readaastringset 関数を用いればアミノ酸配列を読み込めそうだ 31

32 関数の使用法について FASTQ ファイルは読み込めそうだ 32

33 FASTQ 形式ファイル読み込み readdnastringset 関数を用いた読み込み時にオプションを変更することで FASTQ 形式ファイルに対応している 33

34 GC 含量計算の詳細を説明 alphabetfrequency 関数を実行して塩基ごとの出現回数をカウントした結果を hoge に格納している dim 関数は行列 hoge の行数と列数を表示つまり hoge は 7 行 18 列から構成されているということ 34

35 GC 含量計算の詳細を説明 A, C, G, T, および N(A/C/G/T) の出現回数が多いのは当たり前それ以外は M(A/C), R(A/G), W(A/T), S(C/G), といった具合です 35

36 GC 含量計算の詳細を説明入力のシロイヌナズナゲノム配列ファイル (TAIR10_chr_all.fas) 中を検索すると確かに W などの ACGTN 以外の文字が存在します 36

37 GC 含量計算の詳細を説明行列 hoge に対して rowsums 関数を用いて行ごとのカウント数の総和を計算すると染色体ごとの配列長と一致するのは当然です出力 : hoge4.txt 37

38 GC 含量計算の詳細を説明 CG または ACGT のサブセットを抽出してから rowsums 関数を実行 38

39 GC 含量計算の詳細を説明 39

40 パッケージ説明パッケージを個別にインストールする場合使い方の解説記事は PDF のところをクリック例えば 40

41 Biostrings パッケージ中の関数を使いこなせると他の自然言語処理系プログラミング言語 (perl や ruby) を改めて勉強しなくても必要な解析の多くを実行可能 41

42 原著論文引用はお願いします (R で ) 塩基配列解析を利用した証拠もないしアクセスログもとってないので引用や謝辞は必要なし (R で ) 塩基配列解析のことは見なかったことにしても用いた R パッケージや元となるプログラムの原著論文は引用してください by KDT39 42

43 Contents R でゲノム解析シロイヌナズナゲノムの GC 含量計算 multi-fastaファイルの読み込み関数やオプションの利用法パッケージの説明 R でトランスクリプトーム解析シロイヌナズナの RNA-seq データを一通り解析公共 DBからの生データ取得マッピングおよびカウントデータ取得サンプル間クラスタリング発現変動遺伝子 (DEG) 検出 43

44 トランスクリプトーム解析シロイヌナズナの RNA-seq データを一通り R で解析 2 群間比較用 :4 DEX-treated vs. 4 mock-treated 生データ (FASTQファイル) のID:GSE36469 Huang et al., Development, 139: , 2012 生データ取得から発現変動解析までを R のみで実行 44

45 トランスクリプトーム解析シロイヌナズナの RNA-seq データを一通り R で解析 2 群間比較用 :4 DEX-treated vs. 4 mock-treated 生データ (FASTQ ファイル ) の ID:GSE36469 Step1: 生データをダウンロード ( するために必要な ID 情報を取得 ) 45

46 トランスクリプトーム解析シロイヌナズナの RNA-seq データを一通り R で解析 2 群間比較用 :4 DEX-treated vs. 4 mock-treated 生データ (FASTQ ファイル ) の ID:GSE36469 生データをダウンロードするために必要な ID は SRP だということを知る 46

47 トランスクリプトーム解析 SRP を入力として R 上で FASTQ ファイルをダウンロード可能 ( 東大有線 LAN で数時間 ) 実習ではやらないで!!! 47

48 原著論文引用はお願いします CTRLとALTキーを押しながらコードの枠内で左クリックすると全選択できるので積極的に活用 48

49 Step1: 生データのダウンロード中ここでは作業ディレクトリとしてデスクトップ上の SRP を指定している 49

50 Step1: 生データのダウンロード終了後シロイヌナズナの RNA-seq データを一通り R で解析 2 群間比較用 :4 DEX-treated vs. 4 mock-treated ID とサンプル属性 ( ラベル ) との対応関係を知りたい 50

51 トランスクリプトーム解析シロイヌナズナの RNA-seq データを一通り R で解析 2 群間比較用 :4 DEX-treated vs. 4 mock-treated hoge2 実行結果を眺めることで対応付けが可能 51

52 ダウンロード終了後シロイヌナズナの RNA-seq データを一通り R で解析 2 群間比較用 :4 DEX-treated vs. 4 mock-treated ここまでで Step1 生データのダウンロードが完了 52

53 Step2: マッピングおよびカウントデータ取得マッピングに必要な情報 FASTQファイル :8 個の *.fastq.gz リストファイル :srp011435_samplename.txt リファレンスゲノム :TAIR10_chr_all.fas カウントデータ取得に必要な情報遺伝子アノテーションファイル :TAIR10_GFF3_genes.gff 遺伝子ごとにどの染色体のどの座標上に存在するのかなどの情報を含むタブ区切りテキストファイル 53

54 Step2: マッピングおよびカウントデータ取得マッピングに必要な情報リストファイル :srp011435_samplename.txt( 通常はテキストエディタで自作 ) リファレンスゲノム :TAIR10_chr_all.fas(TAIRからダウンロード) カウントデータ取得に必要な情報遺伝子アノテーションファイル :TAIR10_GFF3_genes.gff(TAIRからダウンロード) 必要なファイルを作業ディレクトリに保存 54

55 アノテーションファイル?! 参考 TAIR10 のアノテーションファイル (TAIR10_GFF3_genes.gff) はここからダウンロードしました 55

56 参考 TAIR ウェブインターフェースからアノテーションファイル (TAIR10_GFF3_genes.gff) を取得する際のイメージ 56

57 Step2: マッピングおよびカウントデータ取得 Step2 が二つ存在するがリファレンスとして R パッケージ BSgenome.Athaliana.TAIR.TAIR9 ではなく TAIR10_chr_all.fas を利用するほうで説明します 57

58 最初に description 行の記述を Chr1 や Chr2 に変更した tmp_genome.fasta を作成 58

59 description 行の記述を揃えるのは基本遺伝子アノテーションファイル :TAIR10_GFF3_genes.gff 遺伝子アノテーションファイル中の 1 列目の表記法と同じにするのが基本 59

60 Step2: マッピングおよびカウントデータ取得コード実行後確かに tmp_genome.fasta が作成されている 60

61 Step2: マッピングおよびカウントデータ取得 CTRL と ALT キーを押しながらコードの枠内で左クリックすると全選択できるので積極的に活用 7 時間程度かかるので実行しないで!! 61

62 Step2: マッピングおよびカウントデータ取得私はカウントデータを入力としてその後の各種解析を行います 62

63 Step2: マッピングおよびカウントデータ取得マッピングに必要な情報 FASTQ ファイル :8 個の *.fastq.gz リストファイル :srp011435_samplename.txt リファレンスゲノム :TAIR10_chr_all.fas カウントデータ取得に必要な情報遺伝子アノテーションファイル :TAIR10_GFF3_genes.gff ゲノム上の遺伝子座標情報ファイルを読み込んでいるから遺伝子ごとのカウントデータを取得可能なんですカウントデータファイル :srp011435_count_bowtie_2.txt リストファイル中に記載した任意のサンプル名がカウントデータファイルのヘッダー行となる 63

64 トランスクリプトーム解析シロイヌナズナの RNA-seq データを一通り R で解析 2 群間比較用 :4 DEX-treated vs. 4 mock-treated 生データ (FASTQファイル) のID:GSE36469 Huang et al., Development, 139: , 2012 ここまでで生データ取得からカウントデータ生成まで終了 64

65 トランスクリプトーム解析実験デザイン再確認 2 群間比較用 :4 DEX-treated vs. 4 mock-treated Huang et al., Development, 139: , 2012 群ごとに用いた個体数 (biological replicates) は 2 個体ごとに用いた反復数 (technical replicates) は 2 個体数は 2 群合わせて 4 個体生物アイコン ( 65

66 Step3: サンプル間クラスタリング CTRL と ALT キーを押しながらコードの枠内で左クリックすると全選択できるので積極的に活用 66

67 400 Step3: サンプル間クラスタリング実行結果出力 :srp011435_count_cluster.png

68 Step4: 発現変動遺伝子 (DEG) 同定の前に Technical replicates データのマージを行います 68

69 Step4: 発現変動遺伝子 (DEG) 同定の前に入力 :srp011435_count_bowtie_2.txt 出力 :srp011435_count_bowtie_3.txt Technical replicates データのマージを行います 69

70 Step4: 発現変動遺伝子 (DEG) 同定入力 :srp011435_count_bowtie_3.txt G1 群 G2 群 70

71 390 出力 :srp011435_maplot_bowtie.png

72 出力ファイルの説明正規化後のデータ p-value とその順位 M-A plot の A 値と M 値 q-value FDR 閾値判定結果 q-value < 0.05 を満たす DEG が 1 non-deg が 0 72

73 原著論文引用はお願いします (R で ) 塩基配列解析を利用した証拠もないしアクセスログもとってないので引用や謝辞は必要なし (R で ) 塩基配列解析のことは見なかったことにしても用いた R パッケージや元となるプログラムの原著論文は引用してください by KDT39 73

74 まとめ R でゲノム解析シロイヌナズナゲノムの GC 含量計算 multi-fastaファイルの読み込み関数やオプションの利用法パッケージの説明 R でトランスクリプトーム解析シロイヌナズナの RNA-seq データを一通り解析公共 DBからの生データ取得マッピングおよびカウントデータ取得 Rでいろいろできますサンプル間クラスタリング発現変動遺伝子 (DEG) 検出 74

75 スライド PDF はウェブから取得可能です 75

76 今後の予定 76

77 NGS 速習コース開催東大農 ) 77

シークエンスに基づく比較トランスクリプトーム解析のためのガイドライン構築 ( 代表 ) 新学術領域研究 ( 研究領域提案型 )(H22 年度

78 謝辞共同研究者清水謙多郎先生 ( 東京大学大学院農学生命科学研究科 ) 西山智明先生 ( 金沢大学学際科学実験センター ) 孫建強氏 ( 東京大学大学院農学生命科学研究科大学院生 ) グラント基盤研究 (C)(H24-26 年度 ): シークエンスに基づく比較トランスクリプトーム解析のためのガイドライン構築 ( 代表 ) 新学術領域研究 ( 研究領域提案型 )(H22 年度 -): 非モデル生物におけるゲノム解析法の確立 ( 分担 ; 研究代表者 : 西山智明 ) 挿絵や TCC のロゴなど ( 妻の ) 門田雅世さま作 ( 有能な秘書の ) 三浦文さま作 78

M = log 2 G2 - log 2 G1-2 -1 0 1 2 M-A plot Dudoit et al., Stat.

79 M = log 2 G2 - log 2 G M-A plot Dudoit et al., Stat. Sinica, 12: , 2002 参考 2 群間比較用横軸が全体的な発現レベル縦軸がlog 比からなるプロット名前の由来はおそらく対数の世界での縦軸が引き算 (Minus) 横軸が平均(Average) G1 群 < G2 群 G2 群で高発現 G1 群 = G2 群 G1 群 > G2 群 G1 群で高発現 A = (log 2 G2 + log 2 G1)/2 低発現全体的に高発現 DEG が存在しないデータの M-A plot を眺めることで縦軸の閾値のみに相当する倍率変化を用いた DEG 同定の危険性が分かります 79

80 多重比較問題 :FDR って何? p-value (false positive rate; FPR) 本当は DEG ではないにもかかわらず DEG と判定してしまう確率全遺伝子に占める non-deg の割合 ( 分母は遺伝子総数 ) 例 :10,000 個の non-deg からなる遺伝子を p-value < 0.05 で検定すると 10, = 500 個程度の non-deg を間違って DEG と判定することに相当実際の DEG 検出結果が 900 個だった場合 :500 個は偽物で 400 個は本物と判断実際の DEG 検出結果が 510 個だった場合 :500 個は偽物で 10 個は本物と判断実際の DEG 検出結果が 500 個以下の場合 : 全て偽物と判断 q-value (false discovery rate: FDR) DEG と判定した中に含まれる non-deg の割合 Benjamini and Hochberg J. Roy. Stat. Soc. B, 57: , 1995 DEG 中に占める non-deg の割合 ( 分母は DEG と判定された数 ) non-deg の期待値を計算できれば p 値でも上位 x 個でも DEG と判定する手段はなんでもよい以下は 10,000 遺伝子の検定結果での FDR 計算例 p < を満たす DEG 数が 100 個の場合 :FDR = 10, /100 = 0.1 p < 0.01 を満たす DEG 数が 400 個の場合 :FDR = 10, /400 = 0.25 p < 0.05 を満たす DEG 数が 926 個の場合 :FDR = 10, /926 = 0.54 参考 80

基本的な利用法

(R で ) 塩基配列解析基本的な利用法 Macintosh 版到達目標 : このスライドに書かれている程度のことは自在にできるようにしてエラーへの対処法を身につける 1. 必要なパッケージのインストールが正しくできているかどうかの自力での判定および個別のパッケージのインストール 2. 作業ディレクトリの変更 3. テキストエディタで自在に入出力ファイル名の変更 ( どんなファイル名のものがどこに生成されるかという全体像の把握

ゲノム情報解析基礎 ～ Rで塩基配列解析 ～

ゲノム情報解析基礎～ Rで塩基配列解析～