ネット接続できないヒトも ダブルクリックでローカルに r_seq.html を起動可能です 実習は デスクトップ上にある hoge フォルダの中身が以下の状態を想定して行います (R で ) 塩基配列解析の利用法 : GC 含量計算から発現変動解析まで東京大学 大学院農学生命科学研究科アグリバイオインフォマティクス教育研究プログラム門田幸二 ( かどたこうじ ) kadota@iu.a.u-tokyo.ac.jp http://www.iu.a.u-tokyo.ac.jp/~kadota/ 1
参考ウェブページ (R で ) 塩基配列解析の利用法の紹介 2
前提条件 以下の手順通りに R および必要なパッケージのインストールが完了しているという前提です 3
自己紹介 2002 年 3 月 東京大学 大学院農学生命科学研究科 応用生命工学専攻博士課程修了 学位論文 : cdna マイクロアレイを用いた遺伝子発現解析手法の開発 ( 指導教官 : 清水謙多郎教授 ) 2002/4/1~ 産総研 生命情報科学研究センター (CBRC) 産総研特別研究員 マイクロアレイ解析手法開発 2003/11/1~ 放医研 先端遺伝子発現研究センター研究員 一次元電気泳動波形解析手法開発 2005/2/16~ 東京大学 大学院農学生命科学研究科 アグリバイオインフォマティクスプログラム マイクロアレイ解析手法開発 RNA-seq データ解析手法開発 ( トランスクリプトーム解析周辺の ) 手法開発系のヒトです 4
講義風景 ( 平成 26 年度 ) 5
Contents R でゲノム解析 シロイヌナズナゲノムの GC 含量計算 multi-fastaファイルの読み込み 関数やオプションの利用法 パッケージの説明 R でトランスクリプトーム解析 シロイヌナズナの RNA-seq データを一通り解析 公共 DBからの生データ取得 マッピングおよびカウントデータ取得 サンプル間クラスタリング 発現変動遺伝子 (DEG) 検出 6
植物グローバルなので 例 : シロイヌナズナ (Arabidopsis thaliana) ゲノム配列決定 (chr1-5, chrc, chrm) 1 番染色体 :Theologis et al., Nature, 408: 816-820, 2000 2 番染色体 :Lin et al., Nature, 402: 761-768, 1999 3 番染色体 :Salanoubat et al., Nature, 408: 820-822, 2000 トランスクリプトーム配列 (cdna 配列 ) 決定 アノテーション :Seki et al., Science, 296: 141-145, 2002 まとめサイト The Arabidopsis Information Resource (TAIR) Lamesch et al., Nucleic Acids Res., 40: D1202-1210, 2012 http://www.arabidopsis.org/ R でゲノム解析が可能です 7
アノテーションファイル?! TAIR10 のゲノム配列ファイル (TAIR10_chr_all.fas) はここからダウンロードしました 8
アノテーションファイル?! 参考 TAIR ウェブインターフェースから取得する際のイメージ 9
参考 TAIR10 のゲノム配列ファイル (TAIR10_chr_all.fas) はこんな感じです 10
multi-fasta ファイル?! > のヘッダ行 塩基またはアミノ酸配列 が複数 (multi) 個からなるファイルのこと 11
R で配列長と GC 含量計算 最新のゲノム配列を読み込んで GC 含量を計算することができます TAIR10_chr_all.fas は hoge フォルダ中にあり 12
R の起動 デスクトップにある hoge フォルダ中のファイルを解析 13
作業ディレクトリの変更 Windows(C:) となっている場合もあるが 気にしない この部分はひとそれぞれ 1 2 3 4 5 6 7 14
getwd() と打ち込んで確認 15
基本はコピペ 1 2013 年 7 月以降のリニューアルで コードのコピーがやりずらくなっています CTRL と ALT キーを押しながらコードの枠内で左クリックすると 全選択できます 2 1 一連のコマンド群をコピーして 2R Console 画面上でペースト 16
実行結果 実行前の hoge フォルダ 実行後の hoge フォルダ 出力 : hoge4.txt 17
ありがちなミス 1 作業ディレクトリの変更を忘れているため 入力ファイルの読み込み段階でエラーとなる 18
ありがちなミス 2 必要な入力ファイルが作業ディレクトリ中に存在しない 19
ありがちなミス 3 出力予定のファイル名と同じものを別のプログラムで開いているため最後の write.table 関数のところでエラーが出る 20
ありがちなミス 4 実行スクリプトをコピーする際 最後の行のところで改行を含ませずに R Console 画面上でペーストしたため 最後のコマンドが実行されない ( 出力ファイルが生成されない ) 21
R で配列長と GC 含量計算 出力 : hoge4.txt 原著論文中の数値 ちゃんと似た結果が得られています 22
詳細を説明 出力 : hoge4.txt 入力と出力ファイル名を指定しているところ 23
詳細を説明 GC 含量計算をしたいときには Biostrings というパッケージを読み込む必要があります この作業を行っておかないと 例えば multi-fasta ファイルを読み込むための readdnastringset 関数を利用できません 24
詳細を説明 readdnastringset 関数を用いて 1in_f で指定した入力ファイルを 2fasta 形式で読み込んだ結果を 3fasta というオブジェクト名で格納しています 3 1 2 25
詳細を説明 fasta オブジェクトは確かに multi-fasta ファイル中の情報を適切に読み込めている 26
色についての説明 27
色についての説明 単に確認しているだけなので灰色になっている 28
浮かんでくる疑問 FASTA 形式以外にどんな形式を読み込めるのだろう?FASTQ 形式は読み込めるの? DNA 配列以外にどんな配列を読み込めるのだろう? アミノ酸配列は? 29
関数の使用法について? 関数名 で使用法を記したウェブページが開く ページの下のほうに ( 大抵の場合 ) 使用例が掲載されている 使用法既知の関数のマニュアルをいくつか読んで 慣れておく 30
関数の使用法について FASTQ ファイルは読み込めそうだ readaastringset 関数を用いればアミノ酸配列を読み込めそうだ 31
関数の使用法について FASTQ ファイルは読み込めそうだ 32
FASTQ 形式ファイル読み込み readdnastringset 関数を用いた読み込み時に オプションを変更することで FASTQ 形式ファイルに対応している 33
GC 含量計算の詳細を説明 alphabetfrequency 関数を実行して塩基ごとの出現回数をカウントした結果を hoge に格納している dim 関数は行列 hoge の行数と列数を表示 つまり hoge は 7 行 18 列から構成されているということ 34
GC 含量計算の詳細を説明 A, C, G, T, および N(A/C/G/T) の出現回数が多いのは当たり前 それ以外は M(A/C), R(A/G), W(A/T), S(C/G), といった具合です 35
GC 含量計算の詳細を説明 入力のシロイヌナズナゲノム配列ファイル (TAIR10_chr_all.fas) 中を検索すると 確かに W などの ACGTN 以外の文字が存在します 36
GC 含量計算の詳細を説明 行列 hoge に対して rowsums 関数を用いて行ごとのカウント数の総和を計算すると 染色体ごとの配列長と一致するのは当然です 出力 : hoge4.txt 37
GC 含量計算の詳細を説明 CG または ACGT のサブセットを抽出してから rowsums 関数を実行 38
GC 含量計算の詳細を説明 39
パッケージ説明 パッケージを個別にインストールする場合 使い方の解説記事は PDF のところをクリック 例えば 40
Biostrings パッケージ中の関数を使いこなせると 他の自然言語処理系プログラミング言語 (perl や ruby) を改めて勉強しなくても必要な解析の多くを実行可能 41
原著論文引用はお願いします (R で ) 塩基配列解析 を利用した証拠もないしアクセスログもとってないので引用や謝辞は必要なし (R で ) 塩基配列解析のことは見なかったことにしても 用いた R パッケージや元となるプログラムの原著論文は引用してください by KDT39 42
Contents R でゲノム解析 シロイヌナズナゲノムの GC 含量計算 multi-fastaファイルの読み込み 関数やオプションの利用法 パッケージの説明 R でトランスクリプトーム解析 シロイヌナズナの RNA-seq データを一通り解析 公共 DBからの生データ取得 マッピングおよびカウントデータ取得 サンプル間クラスタリング 発現変動遺伝子 (DEG) 検出 43
トランスクリプトーム解析 シロイヌナズナの RNA-seq データを一通り R で解析 2 群間比較用 :4 DEX-treated vs. 4 mock-treated 生データ (FASTQファイル) のID:GSE36469 Huang et al., Development, 139: 2161-2169, 2012 生データ取得から発現変動解析までを R のみで実行 44
トランスクリプトーム解析 シロイヌナズナの RNA-seq データを一通り R で解析 2 群間比較用 :4 DEX-treated vs. 4 mock-treated 生データ (FASTQ ファイル ) の ID:GSE36469 Step1: 生データをダウンロード ( するために必要な ID 情報を取得 ) 45
トランスクリプトーム解析 シロイヌナズナの RNA-seq データを一通り R で解析 2 群間比較用 :4 DEX-treated vs. 4 mock-treated 生データ (FASTQ ファイル ) の ID:GSE36469 生データをダウンロードするために必要な ID は SRP011435 だということを知る 46
トランスクリプトーム解析 SRP011435 を入力として R 上で FASTQ ファイルをダウンロード可能 ( 東大有線 LAN で数時間 ) 実習ではやらないで!!! 47
原著論文引用はお願いします CTRLとALTキーを押しながらコードの枠内で左クリックすると全選択できるので積極的に活用 48
Step1: 生データのダウンロード中 ここでは作業ディレクトリとして デスクトップ上の SRP011435 を指定している 49
Step1: 生データのダウンロード終了後 シロイヌナズナの RNA-seq データを一通り R で解析 2 群間比較用 :4 DEX-treated vs. 4 mock-treated ID とサンプル属性 ( ラベル ) との対応関係を知りたい 50
トランスクリプトーム解析 シロイヌナズナの RNA-seq データを一通り R で解析 2 群間比較用 :4 DEX-treated vs. 4 mock-treated hoge2 実行結果を眺めることで対応付けが可能 51
ダウンロード終了後 シロイヌナズナの RNA-seq データを一通り R で解析 2 群間比較用 :4 DEX-treated vs. 4 mock-treated ここまでで Step1 生データのダウンロードが完了 52
Step2: マッピングおよびカウントデータ取得 マッピングに必要な情報 FASTQファイル :8 個の *.fastq.gz リストファイル :srp011435_samplename.txt リファレンスゲノム :TAIR10_chr_all.fas カウントデータ取得に必要な情報 遺伝子アノテーションファイル :TAIR10_GFF3_genes.gff 遺伝子ごとに どの染色体のどの座標上に存在するのかなどの情報を含むタブ区切りテキストファイル 53
Step2: マッピングおよびカウントデータ取得 マッピングに必要な情報 リストファイル :srp011435_samplename.txt( 通常はテキストエディタで自作 ) リファレンスゲノム :TAIR10_chr_all.fas(TAIRからダウンロード) カウントデータ取得に必要な情報 遺伝子アノテーションファイル :TAIR10_GFF3_genes.gff(TAIRからダウンロード) 必要なファイルを作業ディレクトリに保存 54
アノテーションファイル?! 参考 TAIR10 のアノテーションファイル (TAIR10_GFF3_genes.gff) はここからダウンロードしました 55
参考 TAIR ウェブインターフェースからアノテーションファイル (TAIR10_GFF3_genes.gff) を取得する際のイメージ 56
Step2: マッピングおよびカウントデータ取得 Step2 が二つ存在するが リファレンスとして R パッケージ BSgenome.Athaliana.TAIR.TAIR9 ではなく TAIR10_chr_all.fas を利用するほうで説明します 57
最初に description 行の記述を Chr1 や Chr2 に変更した tmp_genome.fasta を作成 58
description 行の記述を揃えるのは基本 遺伝子アノテーションファイル :TAIR10_GFF3_genes.gff 遺伝子アノテーションファイル中の 1 列目の表記法と同じにするのが基本 59
Step2: マッピングおよびカウントデータ取得 コード実行後 確かに tmp_genome.fasta が作成されている 60
Step2: マッピングおよびカウントデータ取得 CTRL と ALT キーを押しながらコードの枠内で左クリックすると全選択できるので積極的に活用 7 時間程度かかるので実行しないで!! 61
Step2: マッピングおよびカウントデータ取得 私はカウントデータを入力としてその後の各種解析を行います 62
Step2: マッピングおよびカウントデータ取得 マッピングに必要な情報 FASTQ ファイル :8 個の *.fastq.gz リストファイル :srp011435_samplename.txt リファレンスゲノム :TAIR10_chr_all.fas カウントデータ取得に必要な情報 遺伝子アノテーションファイル :TAIR10_GFF3_genes.gff ゲノム上の遺伝子座標情報ファイルを読み込んでいるから遺伝子ごとのカウントデータを取得可能なんです カウントデータファイル :srp011435_count_bowtie_2.txt リストファイル中に記載した任意のサンプル名がカウントデータファイルのヘッダー行となる 63
トランスクリプトーム解析 シロイヌナズナの RNA-seq データを一通り R で解析 2 群間比較用 :4 DEX-treated vs. 4 mock-treated 生データ (FASTQファイル) のID:GSE36469 Huang et al., Development, 139: 2161-2169, 2012 ここまでで 生データ取得からカウントデータ生成まで終了 64
トランスクリプトーム解析 実験デザイン再確認 2 群間比較用 :4 DEX-treated vs. 4 mock-treated Huang et al., Development, 139: 2161-2169, 2012 群ごとに用いた個体数 (biological replicates) は 2 個体ごとに用いた反復数 (technical replicates) は 2 個体数は 2 群合わせて 4 個体 生物アイコン (http://biosciencedbc.jp/taxonomy_icon/taxonomy_icon.cgi) 65
Step3: サンプル間クラスタリング CTRL と ALT キーを押しながらコードの枠内で左クリックすると全選択できるので積極的に活用 66
400 Step3: サンプル間クラスタリング実行結果 出力 :srp011435_count_cluster.png 500 67
Step4: 発現変動遺伝子 (DEG) 同定の前に Technical replicates データのマージを行います 68
Step4: 発現変動遺伝子 (DEG) 同定の前に 入力 :srp011435_count_bowtie_2.txt 出力 :srp011435_count_bowtie_3.txt Technical replicates データのマージを行います 69
Step4: 発現変動遺伝子 (DEG) 同定 入力 :srp011435_count_bowtie_3.txt G1 群 G2 群 70
390 出力 :srp011435_maplot_bowtie.png 430 71
出力ファイルの説明 正規化後のデータ p-value とその順位 M-A plot の A 値と M 値 q-value FDR 閾値判定結果 q-value < 0.05 を満たす DEG が 1 non-deg が 0 72
原著論文引用はお願いします (R で ) 塩基配列解析 を利用した証拠もないしアクセスログもとってないので引用や謝辞は必要なし (R で ) 塩基配列解析のことは見なかったことにしても 用いた R パッケージや元となるプログラムの原著論文は引用してください by KDT39 73
まとめ R でゲノム解析 シロイヌナズナゲノムの GC 含量計算 multi-fastaファイルの読み込み 関数やオプションの利用法 パッケージの説明 R でトランスクリプトーム解析 シロイヌナズナの RNA-seq データを一通り解析 公共 DBからの生データ取得 マッピングおよびカウントデータ取得 Rでいろいろできます サンプル間クラスタリング 発現変動遺伝子 (DEG) 検出 74
スライド PDF はウェブから取得可能です 75
今後の予定 76
NGS 速習コース開催 (9/1~12@ 東大農 ) 77
謝辞 共同研究者清水謙多郎先生 ( 東京大学 大学院農学生命科学研究科 ) 西山智明先生 ( 金沢大学 学際科学実験センター ) 孫建強氏 ( 東京大学 大学院農学生命科学研究科 大学院生 ) グラント 基盤研究 (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. Sinica, 12: 111-139, 2002 参考 2 群間比較用 横軸が全体的な発現レベル 縦軸がlog 比からなるプロット 名前の由来は おそらく対数の世界での縦軸が引き算 (Minus) 横軸が平均(Average) G1 群 < G2 群 G2 群で高発現 G1 群 = G2 群 G1 群 > G2 群 G1 群で高発現 1 2 3 4 5 A = (log 2 G2 + log 2 G1)/2 低発現 全体的に 高発現 DEG が存在しないデータの M-A plot を眺めることで 縦軸の閾値のみに相当する倍率変化を用いた DEG 同定の危険性が分かります 79
多重比較問題 :FDR って何? p-value (false positive rate; FPR) 本当は DEG ではないにもかかわらず DEG と判定してしまう確率 全遺伝子に占める non-deg の割合 ( 分母は遺伝子総数 ) 例 :10,000 個の non-deg からなる遺伝子を p-value < 0.05 で検定すると 10,000 0.05 = 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: 289-300, 1995 DEG 中に占める non-deg の割合 ( 分母は DEG と判定された数 ) non-deg の期待値を計算できれば p 値でも上位 x 個でも DEG と判定する手段はなんでもよい 以下は 10,000 遺伝子の検定結果での FDR 計算例 p < 0.001 を満たす DEG 数が 100 個の場合 :FDR = 10,000 0.001/100 = 0.1 p < 0.01 を満たす DEG 数が 400 個の場合 :FDR = 10,000 0.01/400 = 0.25 p < 0.05 を満たす DEG 数が 926 個の場合 :FDR = 10,000 0.05/926 = 0.54 参考 80