機能ゲノム学

Size: px

Start display at page:

Download "機能ゲノム学"

こうじのあき
5 years ago
Views:

1 USB メモリ中の hoge フォルダをデスクトップにコピーしておいてください機能ゲノム学第 4 回前回 (5/26) の hoge フォルダがデスクトップに残っているかもしれないのでご注意ください大学院農学生命科学研究科アグリバイオインフォマティクス教育研究プログラム門田幸二 ( かどたこうじ ) kadota@iu.a.u-tokyo.ac.jp 1

2 講義予定第 1 回 (2015 年 5 月 12 日 ) 原理各種データベース生データ取得教科書の 1.2 節 2.2 節周辺第 2 回 (2015 年 5 月 19 日 ) 遺伝子発現行列作成 ( データ正規化 ) クラスタリング ( データ変換や距離の定義など ) 教科書の 3.2 節周辺第 3 回 (2015 年 5 月 26 日 ) 実験デザイン発現変動解析 ( 多重比較問題 ) M-A plot 教科書の 3.2 節と 4.2 節周辺第 4 回 (2015 年 6 月 9 日 ) 機能解析 (Gene Ontology 解析やパスウェイ解析 ) 細胞中で発現している全転写物 ( トランスクリプトーム ) の解析技術はマイクロアレイから次世代シーケンサ ( RNA-seq) に移行しつつありますしかし RNA-seq データ解析の多くはマイクロアレイの知識を前提としています本科目ではマイクロアレイデータを主な例として各種トランスクリプトーム解析手法について解説します教科書 2

3 Contents デザイン行列の意味を理解 ( 教科書 p ) limma パッケージを用いた 2 群間比較のおさらい limma パッケージを用いた 3 群間比較 ( 反復あり ) 反復なし多群間比較 ( 教科書 p ) limma パッケージを用いた 3 群間比較 ( 反復なし ) TCC パッケージ中の ROKU 法を用いた特異的発現遺伝子検出機能解析 ( 遺伝子セット解析 ) 基本的な考え方前処理 MSigDB からの遺伝子セット情報 (gmt 形式ファイル ) 取得 ID 変換 (probe ID gene symbol) GSA パッケージを用いた遺伝子セット解析 3

4 アレイデータ Affymetrix GeneChip Ge et al., Genomics, 86: , 2005 hoge フォルダ中に 3 つの前処理法の実行結果ファイルがあります MAS5 (data_mas_*.txt) RMA (data_rma_*.txt) RMX (data_rob_*.txt) GSE2361 GPL96 (Affymetrix Human Genome U133A Array) 22,283 probesets ヒト 36 サンプル :Heart ( 心臓 ) Thymus ( 胸腺 ) Spleen ( 脾臓 ) Ovary ( 卵巣 ) Kidney ( 腎臓 ) Skeletal Muscle ( 骨格筋 ) Pancreas ( 膵臓 ) Prostate ( 前立腺 ) Nakai et al., Biosci Biotechnol Biochem., 72: , 2008 GSE7623 GPL1355 (Affymetrix Rat Genome Array) 31,099 probesets ラット 24 サンプル :Brown adipose tissue ( 褐色脂肪組織 ; BAT)8 サンプル White adipose tissue ( 白色脂肪組織 ; WAT)8 サンプル Liver ( 肝臓 ; LIV)8 サンプル BAT 8 サンプル : 通常 (BAT_fed) 4 サンプル vs. 24 時間絶食 (BAT_fas) 4 サンプル WAT 8 サンプル : 通常 (WAT_fed) 4 サンプル vs. 24 時間絶食 (WAT_fas) 4 サンプル LIV 8 サンプル : 通常 (LIV_fed) 4 サンプル vs. 24 時間絶食 (LIV_fas) 4 サンプル Kamei et al., PLoS One, 8: e65732, 2013 GSE30533 GPL1355 (Affymetrix Rat Genome Array) 31,099 probesets ラット 10 サンプル : 全て Liver ( 肝臓 ) サンプル iron-deficient diet (Iron_def) 5 サンプル vs. control diet (Control) 5 サンプル 4

5 アレイデータ GSE7623 (Nakai et al., 2008) の対数変換後のデータ data_mas_en.txt data_mas_jp.txt 5

6 2 群間比較 (limma) 教科書 p167- Nakai et al., Biosci Biotechnol Biochem., 72: , 2008 GSE7623 データを用い様々な 2 群間比較を行いクラスタリング結果と DEG 検出結果の関係をみてみよう GSE7623 GPL1355 (Affymetrix Rat Genome Array) 31,099 probesets ラット 24 サンプル :Brown adipose tissue ( 褐色脂肪組織 ; BAT)8 サンプル White adipose tissue ( 白色脂肪組織 ; WAT)8 サンプル Liver ( 肝臓 ; LIV)8 サンプル BAT 8 サンプル : 通常 (BAT_fed) 4 サンプル vs. 24 時間絶食 (BAT_fas) 4 サンプル WAT 8 サンプル : 通常 (WAT_fed) 4 サンプル vs. 24 時間絶食 (WAT_fas) 4 サンプル LIV 8 サンプル : 通常 (LIV_fed) 4 サンプル vs. 24 時間絶食 (LIV_fas) 4 サンプル 1 2 rcode_clustering_png.txt の実行結果 1 肝臓と脂肪間で大きく二つのクラスターに分かれている 2 脂肪の中でも白色脂肪と褐色脂肪に分かれている 3 褐色脂肪は空腹 (24 時間絶食 ) と満腹 ( 通常 ) できれいに分かれている 3 6

7 2 群間比較 (limma) 教科書 p167-4 WAT_fed vs. 4 BAT_fed G1 群 G2 群 7

8 2 群間比較 (limma) 教科書 p167- サブセット抽出のための情報はここで与えている G1 群 G2 群 rcode_limma_4vs4.txt 8

9 2 群間比較 (limma) 教科書 p167- rcode_limma_4vs4.txt 解析したいサブセットに正しくできていることがわかります 9

10 2 群間比較 (limma) 教科書 p167- rcode_limma_4vs4.txt design オブジェクトが ( 実験 ) デザイン行列ですこの行列の 2 列目が G1 群と G2 群がどれに相当するかを表すクラスラベル情報であることもわかります 10

11 2 群間比較 (limma) rcode_limma_4vs4.txt 1limma 実行後の p-value 情報はベクトル形式ではなく行列形式になっていることに注意そしてその列数はデザイン行列 design の列数と同じ 2out$p.value 行列の 2 列目の情報が解析結果に相当

12 2 群間比較 (limma) 1limma 実行後のp-value 情報はベクトル形式ではなく行列形式になっていることに注意そしてその列数はデザイン行列 rcode_limma_4vs4.txt designの列数と同じしつこく示してるだけ 12

13 Contents デザイン行列の意味を理解 ( 教科書 p ) limma パッケージを用いた 2 群間比較のおさらい limma パッケージを用いた 3 群間比較 ( 反復あり ) 反復なし多群間比較 ( 教科書 p ) limma パッケージを用いた 3 群間比較 ( 反復なし ) TCC パッケージ中の ROKU 法を用いた特異的発現遺伝子検出機能解析 ( 遺伝子セット解析 ) 基本的な考え方前処理 MSigDB からの遺伝子セット情報 (gmt 形式ファイル ) 取得 ID 変換 (probe ID gene symbol) GSA パッケージを用いた遺伝子セット解析 13

14 3 群間比較 (limma) 赤枠内では分散分析 (ANOVA) 的な比較するどこかの群間で発現変動している遺伝子の同定について記述しています機能ゲノム学では 1 の教科書に沿った内容で進めます実験デザイン行列だとかコントラスト行列だとか様々な記述形式がありますが 2 応用 1 の例題を一通り行うと感覚がつかめると思います農学生命情報科学特論 I の RNA-seq カウントデータの統計解析のところで説明する予定です

15 3 群間比較 (limma) 教科書 p デザイン行列を理解しておかないと多群間比較時に困ります G1 群 G2 群 G3 群 rcode_limma_4vs4vs4.txt 15

16 3 群間比較 (limma) 教科書 p 解析したいサブセットに正しくできていることがわかります rcode_limma_4vs4vs4.txt 16

17 3 群間比較 (limma) 教科書 p rcode_limma_4vs4vs4.txt 実験デザイン行列には様々が記述の仕方があって難解ですが慣れ以外にはありません 17

18 3 群間比較 (limma) 教科書 p rcode_limma_4vs4vs4.txt levels 関数を用いて合理的にグループラベル情報を抽出しデザイン行列 design の列名情報を変更し取扱いやすくしています 18

3 群間比較 (limma) 教科書 p180-182 rcode_limma_4vs4vs4.

19 3 群間比較 (limma) 教科書 p rcode_limma_4vs4vs4.txt デザイン行列の列名を変更して取扱いやすくしておかないとこの部分での指定時にややこしいことになるここでは 3 種類の 2 群間比較を行うようにしているここで作成しているのはコントラスト行列というもの比較したいものを作成した行列という位置づけ 19

20 3 群間比較 (limma) 教科書 p rcode_limma_4vs4vs4.txt 3 種類の 2 群間比較を行うようにしたコントラスト行列 contrast を入力としているので DEG 検出結果として 31,099 行 3 列からなる p- value 行列が得られることになる 20

21 3 群間比較 (limma) 教科書 p rcode_limma_4vs4vs4.txt apply 関数を用いて列ごと (MARGIN=2) に q-value を計算している 21

3 群間比較 (limma) 教科書 p180-182 FDR 1% という閾値を満たす遺伝子数は G1 vs.

22 3 群間比較 (limma) 教科書 p FDR 1% という閾値を満たす遺伝子数は G1 vs. G2 が 2,418 個 G1 vs. G3 が 8,119 個そして G2 vs. G3 が 7,471 個という結果 22

23 3 群間比較 (limma) 教科書 p G1vsG2 の DEG 数が他に比べて少ないので妥当 G1 群 G2 群 G3 群 23

24 Contents デザイン行列の意味を理解 ( 教科書 p ) limma パッケージを用いた 2 群間比較のおさらい limma パッケージを用いた 3 群間比較 ( 反復あり ) 反復なし多群間比較 ( 教科書 p ) limma パッケージを用いた 3 群間比較 ( 反復なし ) TCC パッケージ中の ROKU 法を用いた特異的発現遺伝子検出機能解析 ( 遺伝子セット解析 ) 基本的な考え方前処理 MSigDB からの遺伝子セット情報 (gmt 形式ファイル ) 取得 ID 変換 (probe ID gene symbol) GSA パッケージを用いた遺伝子セット解析 24

25 教科書反復なし 3 群間比較 (limma) (biological) replicates がないデータの場合 limma は基本的に適用不可 G1 群 G2 群 G3 群 rcode_limma_1vs1vs1.txt 25

26 教科書反復なし 3 群間比較 (limma) (biological) replicates がないデータの場合 limma は基本的に適用不可 26

27 反復なし多群間比較 (TCC) 入力教科書 ROKU 入出力のイメージ出力は入力の対応する位置に特異的組織でない部分には 0 特異的高発現が 1 特異的低発現が -1 と返す 1 出力の右側は全体的な組織特異度の高さでランキングした結果を返す黒枠内の遺伝子 (gene2, 7, 8) は組織特異的遺伝子ではない 1 出力 27

28 反復なし多群間比較 (TCC) 教科書高発現は 1 低発現は -1 入力出力 28

29 教科書反復なし多群間比較 (TCC) 入出力のイメージはこの例題実行結果 modh 列の値がデータ変換後のエントロピー値に相当 ( の講義でほんの少しだけ触れました ) 29

30 IC Sequence logos position iの情報量 IC ( N) H( x ) i log 2 2 参考 i Sequence logos はあるポジションに特定の塩基が濃縮されている状態をうまく表すためにエントロピーを内部的に計算しているだけです水色の枠内がエントロピー p 1,4 = 90% p 5,3 = 50% p 5,1 = 50% Schneider and Stephens, Nucleic Acids Res., 18(20): ,

31 エントロピー遺伝子 i のエントロピー H(x ) i 参考組織特異的発現遺伝子検出にエントロピーを最初に応用したのは Schug ら (2005) N p j ij log pij pij x 1 2 ( ), where ij / N j 1 x ij 組織特異的遺伝子は低いエントロピー N: 組織数 (jの数) = 8 Hの取りうる範囲 :0 H log 2 N 0 H 3 そうでないものは高い値 Schug et al., Genome Biol., 6: R33,

教科書 4.2.3 ROKU を実行 Affymetrix GeneChip Ge et al.

GPL96 (Affymetrix Human Genome U133A Array) 22,283 probesets ヒト 36 サンプル :Heart (

32 教科書 ROKU を実行 Affymetrix GeneChip Ge et al., Genomics, 86: , 2005 hoge GSE2361 フォルダ中のデータを用いて ROKU を実行してみよう GSE2361 GPL96 (Affymetrix Human Genome U133A Array) 22,283 probesets ヒト 36 サンプル :Heart ( 心臓 ) Thymus ( 胸腺 ) Spleen ( 脾臓 ) Ovary ( 卵巣 ) Kidney ( 腎臓 ) Skeletal Muscle ( 骨格筋 ) Pancreas ( 膵臓 ) Prostate ( 前立腺 ) 32

Additional file 1 と同じです GSE2361 GPL96 (Affymetrix Human Genome U133A Array) 22,283 probesets

33 教科書 ROKU を実行 Affymetrix GeneChip Ge et al., Genomics, 86: , 2005 赤枠の RMA 前処理後のデータを入力として ROKU を実行した結果 ( の一部 ) が基本的に ROKU 論文の Additional file 1 と同じです GSE2361 GPL96 (Affymetrix Human Genome U133A Array) 22,283 probesets ヒト 36 サンプル :Heart ( 心臓 ) Thymus ( 胸腺 ) Spleen ( 脾臓 ) Ovary ( 卵巣 ) Kidney ( 腎臓 ) Skeletal Muscle ( 骨格筋 ) Pancreas ( 膵臓 ) Prostate ( 前立腺 ) 33

34 ROKU を実行教科書テンプレートスクリプトをテキストエディタにコピペして入力ファイル名部分を変更し R Console 画面上でコピペ約 2 分かかります 34

35 ROKU を実行教科書 ROKU 論文の Additional file 1 と若干結果が違いますが R やパッケージのバージョンの違い ROKU 内部で用いている階乗計算部分が原著論文では近似式 (Staring) でしたが TCC に移植する際に別の関数に置換したことなどが原因と考えられます 35

36 Tips ROKU の詳細について記述していますまたリクエストのあった 1WAD 法の数式や 2 前処理法と DEG 検出法の組合せのガイドラインについてもこの PDF に記載あり教科書 p170 にも記載あり 36

37 Contents デザイン行列の意味を理解 ( 教科書 p ) limma パッケージを用いた 2 群間比較のおさらい limma パッケージを用いた 3 群間比較 ( 反復あり ) 反復なし多群間比較 ( 教科書 p ) limma パッケージを用いた 3 群間比較 ( 反復なし ) TCC パッケージ中の ROKU 法を用いた特異的発現遺伝子検出機能解析 ( 遺伝子セット解析 ) 基本的な考え方前処理 MSigDB からの遺伝子セット情報 (gmt 形式ファイル ) 取得 ID 変換 (probe ID gene symbol) GSA パッケージを用いた遺伝子セット解析 37

38 機能解析 Gene Ontology (GO) 解析 ( 発現に差のある GO term を探索 ) 基本 3 カテゴリ (Cellular Component (CC), Molecular Function (MF), Biological Process (BP)) のどれでも可能例 : 肝臓の空腹状態 vs. 満腹状態の GO(BP) 解析の結果脂肪酸 β 酸化関連 GO term (GO: ) が動いていることが分かったパスウェイ解析 ( 発現に差のあるパスウェイを探索 ) KEGG Pathway, BioCarta, Reactome pathway database のどれでも可能例 : 酸化的リン酸化パスウェイ関連遺伝子セットが糖尿病患者で動いていたモチーフ解析 ( 発現に差のあるモチーフを探索 ) 同じ 3 -UTR microrna 結合モチーフをもつ遺伝子セット同じ転写因子結合領域 (TATA-box など ) をもつ遺伝子セット例 :TATA-box をもつ遺伝子セットが G1 群 vs. G2 群比較で動いていた機能解析の実体は遺伝子セットの発現変動解析発現に差のある遺伝子セットを探したいということ 38

39 N =10,000 genes 機能解析 ( 遺伝子セット解析 ) 発現変動遺伝子セット解析手法 (2 群間比較用がほとんど ) N =10,000 個の遺伝子からなる2 群間比較用データこの中に XXX 関連遺伝子がn 個含まれている例 : 酸化的リン酸化 (=XXX) 関連遺伝子が 7(=n) 個含まれている酸化的リン酸化関連遺伝子セットが変動しているかどうかを調べたいという問題を考える G1 群 G2 群 7 個の酸化的リン酸化関連遺伝子の位置 39

40 N =10,000 genes 機能解析 ( 遺伝子セット解析 ) 遺伝子ごとの発現変動の度合いを数値化例 :t- 統計量 log2(g2/g1) 相関係数基本はデフォルトでやるようですが様々な選択肢があります G1 群 G2 群 G1 群 G2 群発現変動遺伝子 (G1 群 > G2 群 ) 変動してない遺伝子発現変動遺伝子 (G1 群 < G2 群 ) 40

41 N =10,000 genes 機能解析 ( 遺伝子セット解析 ) 発現変動順にソート後の酸化的リン酸化関連遺伝子セットのステレオタイプな分布どうやって偏りを評価するのか? 変動している変動してない G1 群 G2 群 G1 群 G2 群 41

42 N =10,000 genes 機能解析 ( 遺伝子セット解析 ) Over-Representation Analysis (ORA) 何らかの手段で決めた上位 X(=1500) 個のうち x 個が酸化的リン酸化関連遺伝子であった基本的な考え方は全遺伝子と上位のサブセットのみで調べたい遺伝子セットの割合が不変という帰無仮説のもとで検定酸化的リン酸化関連遺伝子セット (n =7) が変動していない場合 : x/n X/N (= 1500/10000) 酸化的リン酸化関連遺伝子セット (n =7) が変動している場合 : x/n >> X/N (= 15%) G1 群 G2 群 G1 群 G2 群

43 N =10,000 genes 機能解析 ( 遺伝子セット解析 ) Over-Representation Analysis (ORA) 何らかの手段で決めた上位 X(=1500) 個のうち x 個が酸化的リン酸化関連遺伝子であった 2 2 分割表 (contingency table) に基づく方法超幾何検定やカイ二乗検定が利用されます G1 群 G2 群 G1 群 G2 群 6 XXX= 酸化的リン酸化関連遺伝子セット 43

44 rcode_ora_basic.txt 機能解析 ( 超幾何検定 ) DEG として 1500 個抽出したとき酸化的リン酸化関連遺伝子が 6 個以上含まれる確率として算出 N=10000 個の遺伝子発現データ中に XXX= 酸化的リン酸化関連遺伝子は n=7 個含まれていた上位 X=1500 個の発現変動遺伝子 (DEG) の中に x=6 個の酸化的リン酸化関連遺伝子が含まれていた帰無仮説 : 酸化的リン酸化関連遺伝子の割合は DEG と non-deg 間で差がない 44

45 rcode_ora_basic.txt 機能解析 ( 超幾何検定 )?dhyper マニュアル中の一般的な説明に置き換えるとこんな感じです m=7 個の白いボールと n=9993 個の黒いボールが入った箱があります ( トータルで N=m+n=10,000 個 ) この中から k=1500 個ランダムに取り出したときに x=6 個以上白いボールが含まれる確率を計算しなさい 45

46 rcode_ora_basic.txt 機能解析 ( カイ二乗検定 ) DEG として 1500 個抽出したとき酸化的リン酸化関連遺伝子が 6 個以上含まれる確率として算出 46

47 N =10,000 genes 機能解析 ( 遺伝子セット解析 ) Over-Representation Analysis (ORA) 上位 1500 個のうち酸化的リン酸化関連遺伝子が7 個中 4つ以上含まれていればp < 0.05で検出可能ということを意味する G1 群 G2 群 G1 群 G2 群 p < 0.05 を灰色で示した 47

48 機能解析 ( 遺伝子セット解析 ) Over-Representation Analysis (ORA) GenMAPP (Dahlquist et al., Nature Genet., 31: 19-20, 2002) FatiGO (Al-Shahrour et al., Bioinformatics, 20: , 2004) GOstat (Beissbarth et al., Bioinformatics, 20: , 2004) GOFFA (Sun et al., BMC Bioinformatics, 7 Suppl 2: S23, 2006) agrigo (Du et al., Nucleic Acids Res., 38: W64-W70, 2010) GenMAPP は比較的有名だと思います 48

49 N =10,000 genes 第 1 世代 (ORA) の短所 1 全体的には動いているものの個々の発現変動の度合いが弱い場合に検出困難 2 上位 X 個のX 次第で結果が変わる 3 情報量低下 ( 発現変動の度合いカウント情報 ) もちろん分割表ベースの方法 (ORA) ではない第 2 世代以降の方法があります代表例は Gene Set Enrichment Analysis (GSEA) 3 G1 群 G2 群 2 G1 群 G2 群

50 第 2 世代 (FCS) もちろん分割表ベースの方法 (ORA) ではない第 2 世代以降の方法があります代表例は Gene Set Enrichment Analysis (GSEA) Functional Class Scoring (FCS) 1. 遺伝子ごとの統計量を算出 ( 発現変動の度合いを数値化 ) 例 :t- 統計量 log(g2/g1) 相関係数 2. 目的の遺伝子セット XXX(= 酸化的リン酸化関連遺伝子 ) の偏りを何らかの方法で評価 t 検定 (XXX 中の遺伝子群の統計量 vs. それ以外の遺伝子群の統計量 ) Wilcoxon rank sum test (XXX 中の遺伝子群の発現変動の順位 vs. それ以外 ) XXX 中の n 個の遺伝子群の何らかの要約統計量 S XXX を計算しておき N 個の全遺伝子の中からランダムに n 個を抽出して同じ統計量を計算する ( 例えば 10 万回 ) 10 万回のうち S XXX 以上 ( 大きければ大きいほど発現変動していることを意味する場合 ; その逆のときは以下 ) だった回数 ( 例えば j 回 ) に基づいて p 値 (= j / 100,000) を算出 ( いわゆる gene set permutation というアプローチ ) 本来の G1 群 vs. G2 群のラベル情報を用いて得られた XXX 中の n 個の遺伝子群の何らかの要約統計量 S XXX を計算しておくランダムにラベル情報を入れ替えて同じ統計量を計算することを何回も繰り返して p 値を算出 ( いわゆる Phenotype permutation というアプローチ ) 50

51 N =10,000 genes 第 2 世代 (FCS) 第一世代 (ORA) の欠点が改善 1 全体的には動いているものの個々の発現変動の度合いが弱い場合に検出困難 2 上位 X 個のX 次第で結果が変わる 3 情報量低下 ( 発現変動の度合いカウント情報 ) 遺伝子ごとの log 比で考えると遺伝子を等価に取り扱うのではなく log 比そのものを足し込むことで発現変動の大きなものと小さなものを考慮するようなイメージ 3 G1 群 G2 群 2 G1 群 G2 群

52 第 2 世代 (FCS) Functional Class Scoring (FCS) GSEA (Subramanian et al., PNAS, 102: , 2005) PAGE (Kim and Volsky, BMC Bioinformatics, 6: 144, 2005) sigpathway (Tian et al., PNAS, 102: , 2005) GSA (Efron and Tibshirani, Ann. Appl. Stat., 1: , 2007) GeneTrail (Backes et al., Nucleic Acids Res., 35: W186-W192, 2007) SAM-GS (Dinu et al., BMC Bioinformatics, 8: 242, 2007) 最も有名なのは GSEA です 52

53 Khatri et al., PLoS Comput. Biol., 8(2): e , 2012 突っ込みどころは満載だがそんなことをいってもしょうがない遺伝子セット解析の課題 ( 知識ベースの解析法なので ) 解析対象がアノテーションの情報の豊富な生物種に限定それ以外の生物種はまずは地道にアノテーション情報を増やしていくことが先決 ( ではないだろうか ) アノテーションの解像度を上げる努力も大事アノテーション情報の信頼度が高いとはいえないなんらかの GO term がついていたとしてもその大部分の evidence code が自動でつけられたもの (IEA, inferrred from electronic annotations) である遺伝子セット間の独立性の問題数百個程度の遺伝子セットの中から比較するサンプル間で動いている遺伝子セットはどれか? という解析を遺伝子セット間の独立性を仮定して調べるがそもそも独立ではない (GO term 間の親子関係などから明らか ) いくつくらいの遺伝子セットが動いているのか? という問いに答えるすべがない評価に用いられるよく研究されているデータセットは答えが完全に分かっているものではない (the actual biology is never fully known!) 感度が高いと謳っているだけの方法は ( 全部の遺伝子セットが動いている感度 100%) 53

54 遺伝子セット解析おさらい Gene Ontology (GO) 解析 ( 発現に差のある GO term を探索 ) 基本 3 カテゴリ (Cellular component (CC), Molecular Function (MF), Biological Process (BP)) のどれでも可能例 : 肝臓の空腹状態 vs. 満腹状態の GO(BP) 解析の結果脂肪酸 β 酸化関連 GO term (GO: ) が動いていることが分かったパスウェイ解析 ( 発現に差のあるパスウェイを探索 ) KEGG, BioCarta, Reactome pathway database のどれでも可能例 : 酸化的リン酸化パスウェイ関連遺伝子セットが糖尿病患者で動いていたモチーフ解析 ( 発現に差のあるモチーフを探索 ) Subramanian et al., PNAS, 102: , 2005 同じ 3 -UTR microrna 結合モチーフをもつ遺伝子セット同じ転写因子結合領域 (TATA-box など ) をもつ遺伝子セット例 :TATA-box をもつ遺伝子セットが G1 群対 G2 群比較で動いていたどの遺伝子セットにどの遺伝子が所属しているかという gmt 形式ファイルの取得が第一歩 54

55 Contents デザイン行列の意味を理解 ( 教科書 p ) limma パッケージを用いた 2 群間比較のおさらい limma パッケージを用いた 3 群間比較 ( 反復あり ) 反復なし多群間比較 ( 教科書 p ) limma パッケージを用いた 3 群間比較 ( 反復なし ) TCC パッケージ中の ROKU 法を用いた特異的発現遺伝子検出機能解析 ( 遺伝子セット解析 ) 基本的な考え方前処理 MSigDB からの遺伝子セット情報 (gmt 形式ファイル ) 取得 ID 変換 (probe ID gene symbol) GSA パッケージを用いた遺伝子セット解析 55

56 Subramanian et al., PNAS, 102: , 2005 MSigDB ver. 5.0 H: hallmark gene sets (50 gene sets) c1: positional gene sets (326 gene sets) ヒト染色体の位置ごとの遺伝子セットリストファイル (326 gene sets) c2: curated gene sets (4,725 gene sets) CGP: chemical and genetic perturbations (3,395 gene sets) CP: canonical pathways (1,330 gene sets) CP:BIOCARTA: BioCarta gene sets (217 gene sets) CP:KEGG: KEGG gene sets (186 gene sets) CP:REACTOME: Reactome gene sets (674 gene sets) c3: motif gene sets (836 gene sets) MIR: microrna targets (221 gene sets) TFT: transcription factor targets (615 gene sets) c4: computational gene sets (858 gene sets) CGM: cancer gene neighborhoods (427 gene sets) CM: cancer modules (431 gene sets) c5: gene ontology (GO) gene sets (1,454 gene sets) BP: biological process (825 gene sets) CC: cellular component (233 gene sets) MF: molecular function (396 gene sets) c6: oncogenic signatures gene sets (189 gene sets) c7: immunologic signatures gene sets (1,910 gene sets) MSigDB は Molecular Signature Database の略様々な遺伝子セット解析を行うための gmt 形式ファイルをダウンロード可能です発現変動と関連する KEGG パスウェイを調べたいとき発現変動と関連する BP 中の GO terms を調べたいとき 56

57 Subramanian et al., PNAS, 102: , 2005 MSigDB ver 年 4 月に ver. 4.0 から 5.0 になったようです劇的な違いはないようです c1: positional gene sets (326 gene sets) ヒト染色体の位置ごとの遺伝子セットリストファイル (326 gene sets) c2: curated gene sets (4,722 gene sets) CGP: chemical and genetic perturbations (3,402 gene sets) CP: canonical pathways (1,320 gene sets) CP:BIOCARTA: BioCarta gene sets (217 gene sets) CP:KEGG: KEGG gene sets (186 gene sets) CP:REACTOME: Reactome gene sets (674 gene sets) c3: motif gene sets (836 gene sets) MIR: microrna targets (221 gene sets) TFT: transcription factor targets (615 gene sets) c4: computational gene sets (858 gene sets) CGM: cancer gene neighborhoods (427 gene sets) CM: cancer modules (431 gene sets) c5: gene ontology (GO) gene sets (1,454 gene sets) BP: biological process (825 gene sets) CC: cellular component (233 gene sets) MF: molecular function (396 gene sets) c6: oncogenic signatures gene sets (189 gene sets) c7: immunologic signatures gene sets (1,910 gene sets) 発現変動と関連する KEGG パスウェイを調べたいとき発現変動と関連する BP 中の GO terms を調べたいとき 57

58 MSigDB Subramanian et al., PNAS, 102: , 2005 遺伝子セット解析を行うための gmt 形式ファイルのダウンロード方法はこちら

59 MSigDB Subramanian et al., PNAS, 102: , 発現変動と関連する biological processes (BP) 中の GO terms を調べたいときは 3 黒枠内のいずれかの gmt ファイルを利用

60 gmt ファイル基本的にどれを使っても自由だが利用する R パッケージがどの入力形式を受け付けるかにも依存する経験上 gene symbols を使っておけば間違いないので門田は *.symbols.gmt をいつも利用しています 1 列目 : 遺伝子セット名 2 列目 :URL 3 列目以降 :gene ID or symbol 60

61 GO 解析 GSE7623 (Nakai et al., 2008) の対数変換後のデータを入力として BAT_fed vs. BAT_fas の遺伝子セット解析をやってみよう G1 群 G2 群 61

62 GO 解析 ( 前処理 ) プローブ ID と gene symbol の対応付けを行い同じ gene symbol に複数のプローブ ID が割り当てられる場合は平均値を採用するなどして non-redundant にする ( 折り畳む ; つぶす ;collapse) 作業が必要 G1 群 G2 群 62

63 Contents デザイン行列の意味を理解 ( 教科書 p ) limma パッケージを用いた 2 群間比較のおさらい limma パッケージを用いた 3 群間比較 ( 反復あり ) 反復なし多群間比較 ( 教科書 p ) limma パッケージを用いた 3 群間比較 ( 反復なし ) TCC パッケージ中の ROKU 法を用いた特異的発現遺伝子検出機能解析 ( 遺伝子セット解析 ) 基本的な考え方前処理 MSigDB からの遺伝子セット情報 (gmt 形式ファイル ) 取得 ID 変換 (probe ID gene symbol) GSA パッケージを用いた遺伝子セット解析 63

64 ID 変換 1 教科書 p70-71 遺伝子発現データは公共 DB の GEO から GSE7623 という ID で取得したものだったここからプローブ ID と gene symbol の対応付けを行うためのアノテーションファイルを取得可能 2 64

65 ID 変換教科書 p70-71 プローブ ID と gene symbol からなるアノテーションファイルを取得できています確認時は 2 分程度で終わりましたが hoge フォルダに hoge3_gpl1355.txt を一応置いてあります 65

66 ID 変換エクセルで開くときには注意が必要!1 行 1 列目のところが ID から始まる文字列の場合にこのような現象が起こるようですが基本無視で構いません

67 ID 変換参考編集して保存したい場合にはドラッグ & ドロップで開いてはだめですファイル - 開くでファイルを指定して開くべし! そのまま開くと例えば March2 という gene symbol が日付と認識されてしまうためこれを防ぐ必要があります! 67

68 ID 変換ここではファイルの中身を眺めるだけなので再度ドラッグ & ドロップ 1 回目は失敗しても 2 回目は普通に開けます 68

69 ID 変換 Gene Symbol 列でソートしてみると hoge3_gpl1355.txt data_mas_en.txt 69

70 ID 変換同じ gene symbol を持つプローブ ID が複数存在することがわかる Gene Symbol 列でソート 70

71 ID 変換入力 1:hoge3_GPL1355.txt 入力 2:data_mas_EN.txt マイクロアレイごとに搭載されている遺伝子の種類や重複度が異なるためこの作業は重要出力 :data_mas_en_symbol.txt 71

72 ID 変換 2 つの入力ファイル ( 発現データと変換表 ) から 1 つの出力ファイルが得られます 72

73 ID 変換 hoge フォルダ中の data_mas_en_symbol.txt はこのコードのコピペで作成しています作業ディレクトリに入力ファイルがあることを確認してから実行しましょう rcode_id_conversion.txt 73

74 ID 変換 rcode_id_conversion.txt hoge フォルダ中の data_mas_en_symbol.txt は 14,132 個のユニークな gene symbols からなることがわかります 74

75 Tips:as.matrix プログラムの組み方で速度が結構違います ( データフレーム形式より行列形式のほうが早いらしい ) 孫建強氏作は 1 分門田作は 2 分 ( 爆 ) 75

76 Contents デザイン行列の意味を理解 ( 教科書 p ) limma パッケージを用いた 2 群間比較のおさらい limma パッケージを用いた 3 群間比較 ( 反復あり ) 反復なし多群間比較 ( 教科書 p ) limma パッケージを用いた 3 群間比較 ( 反復なし ) TCC パッケージ中の ROKU 法を用いた特異的発現遺伝子検出機能解析 ( 遺伝子セット解析 ) 基本的な考え方前処理 MSigDB からの遺伝子セット情報 (gmt 形式ファイル ) 取得 ID 変換 (probe ID gene symbol) GSA パッケージを用いた遺伝子セット解析 76

77 GO 解析の準備完了褐色脂肪満腹 vs. 空腹の発現変動に関連した GO Biological Process 遺伝子セットを GSA 法で解析するための前処理が完了入力 1:data_mas_EN_symbol.txt G1 群 G2 群入力 2:c5.bp.v5.0.symbols.gmt 77

78 GSA で GO 解析 data_mas_en_symbol.txt を入力として BAT_fed vs. BAT_fas の GO 解析をやってみよう

79 GSA で GO 解析 FDR 10% の閾値を満たす有意な遺伝子セット数は G1 群で高発現のものが 24 個 G2 群で高発現のものが 4 個だったことがわかるヒトによって若干結果が異なります 79

80 GSA で GO 解析私は結果の評価はできません 80

機能ゲノム学（第6回）

機能ゲノム学（第6回）講義室後ろにある USB メモリ中の hoge フォルダをデスクトップにコピーしておいてください機能ゲノム学第 4 回東京大学大学院農学生命科学研究科アグリバイオインフォマティクス教育研究ユニット門田幸二 kadota@iu.a.u-tokyo.ac.jp 前回 (5/28) のhogeフォルダがデスクトップに残っているかもしれないのでご注意ください Jun 04, 2014 1 講義予定