機能ゲノム学（第6回）

Size: px

Start display at page:

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

あきみしどり
5 years ago
Views:

1 マイクロアレイデータ解析結果の正しい?! 解釈について東京大学大学院農学生命科学研究科アグリバイオインフォマティクス教育研究ユニット門田幸二 ( かどたこうじ ) kadota@iu.a.u-tokyo.ac.jp 1

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

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

3 講義内容アレイデータの正規化 ( 前処理 ) 生データ遺伝子発現行列クラスタリング発現変動遺伝子 (DEG) の同定二群間比較評価基準評価法および (Affymetrix チップの ) ガイドライン多サンプル間比較組織特異的遺伝子機能解析 (GSEA 解析 ) Gene Ontology 解析パスウェイ解析 3

4 4

5 アレイデータの正規化 ( 前処理 ) 実験によって得られた生のシグナル強度をそのまま利用することは普通はやりません二色法 : 蛍光色素 (Cy3 and Cy5) の取り込み効率補正一色法 : シグナルゲイン?! の補正こうであるべき! という仮定を置いてそれを満たすような正規化を行った後のデータを利用する 5

6 プローブセット Affymetrix 製チップ解析戦略遺伝子 i の発現量 S i を n i (n i =11~20) 種類のプローブペアのシグナル強度をもとに計算 5 3 PM PM PM PM PM PM PM PM PM PM PM i,1 i,2 i,3 i,4 i,5 i,6 i,7 i,8 i,9 i,10 i,11, MM, MM, MM, MM, MM, MM, MM, MM, MM i,1, MM, MM i,2 i,3 i,4 i,5 i,6 i,7 i,8 i,9 i,10 i,11 5 CAGAATCATTAGACTATCCGATAAGGAGTACAATCTGA 3 CATTAGACTATCCGATAAGGAGTAC Perfect match (PM i,j ) プローブペア CATTAGACTATCGGATAAGGAGTAC Mismatch (MM i,j ) 25 mer 遺伝子 iの発現量 S i ( summary score or expression index ) 発現量 S i を算出するための様々な前処理法が存在 6

7 Affymetrix 製チップ解析戦略 ( 様々な前処理法 ) MBEI (Li and Wong, PNAS, 98, 31-36, 2001) MAS5 (Hubbell et al., Bioinformatics, 18, , 2002) RMA (Irizarry et al., Biostatistics, 4, , 2003) GCRMA (Wu et al., Tech. Rep., John Hopkins Univ., 2003) PDNN (Zhang et al., Nat. Biotechnol., 21, , 2003) PLIER (Affymetrix, 2004) SuperNorm (Konishi, T., BMC Bioinformatics, 5, 5, 2004) multi-mgmos (Liu et al., Bioinformatics, 21, , 2005) GLA (Zhou and Rocke, Bioinformatics, 21, , 2005) FARMS (Hochreiter et al., Bioinformatics, 22, , 2006) DFW (Chen et al., Bioinformatics, 23, , 2007) Hook (Binder et al., AMB, 3, 11, 2008) 生データ ( PM i, j, MM i, ) j in.cel files バックグラウンド補正 (within-array) 正規化 (crossarray) PM 値の補正 Summarization 発現量 S i 7

8 グローバル正規化仮定 : 各サンプルから測定された mrna の全体量は一定チップ上の遺伝子数が尐ない場合は非現実的だが数千 ~ 数万種類の遺伝子が搭載されているので妥当 ( だろう ) nomalization 2008/7/16 8

9 Quantile 正規化仮定 : 順位が同じならシグナル強度も同じ正規化前正規化後列ごとにソート行ごとの平均を算出対応する行の要素の元の位置に平均値を代入データセット中のサンプル数が変わると結果が変わる 9

10 正規化遺伝子発現行列二群間比較様々な組織 ( 条件 ) 時系列データ A A B B x 1, 1 x 1, 2 x 1, 2 x 1, 2 A A B B x 2, 1 x 2, 2 x 2, 2 x 2, 2 x 1,1 x 2,1 x 1,2 x 2,2 x 1,3 x 2,3 x 1,4 x 2,4 x 1,1 x 2,1 x 1,2 x 2,2 x 1,3 x 2,3 x 1,4 x 2,4 A xi, 1 A B B xi, 2 xi, 2 xi, 2 x i,1 x i,2 x i,3 x i,4 x i,1 x i,2 x i,3 x i,4 A A B B xn, 1 xn, 2 xn, 2 xn, 2 x n,1 x n,2 x n,3 x n,4 x n,1 x n,2 x n,3 x n,4 様々な解析が可能な状態 10

11 手順 1: 前処理法の適用 Affymetrix GeneChip の場合様々な前処理法を適用し複数の遺伝子発現行列データを得る例 )MAS, RMA, qfarms or DFW その他のメーカーの場合メーカー推奨のやり方に従って遺伝子発現行列データ ( 基本的に一つのみ ) を得る理由 : どの前処理法を使うかでサンプル間クラスタリング ( 後述 ) の結果が大きく異なりうるから 11

12 クラスタリングサンプルの属性情報 ( 癌 or 正常など ) を使わずに発現情報のみを用いて発現パターンの類似した遺伝子 ( またはサンプル ) をクラスター ( 群 ) にしていく手法 (Unsupervised learning 二群間比較多サンプル時系列解析 A A B B x 1, 1 x 1, 2 x 1, 2 x 1, 2 A A B B x 2, 1 x 2, 2 x 2, 2 x 2, 2 x 1,1 x 2,1 x 1,2 x 2,2 x 1,3 x 2,3 x 1,4 x 2,4 x 1,1 x 2,1 x 1,2 x 2,2 x 1,3 x 2,3 x 1,4 x 2,4 A xi, 1 A B B xi, 2 xi, 2 xi, 2 x i,1 x i,2 x i,3 x i,4 x i,1 x i,2 x i,3 x i,4 A A B B 2009/08/19 xn, 1 xn, 基礎生物学研究所 2 xn, 2 xn, 2 x n,1 x n,2 x n,3 x n,4 x n,1 x n,2 x n,3 x n,

13 Bittner et al., Nature, 2000 サンプル間クラスタリングの例メラノーマサンプル悪性度の高い癌のサブタイプを発見 13

14 クラスタリング階層的クラスタリング発現パターンの類似した遺伝子を集めて系統樹を作成非階層的クラスタリング K-means クラスタリング K 個のクラスターに分割 (K の数は主観的に決定 ) すると予め指定し各クラスター内の遺伝子 ( サンプル ) 間の距離の総和が最小になるような K 個のクラスターを作成自己組織化マップ (SOM) 主成分分析 (PCA) 14

15 距離 ( 類似度 ) の定義遺伝子 (or サンプル )x と y の発現パターンの距離 D 1) 1 ( ) ( 1 1 ) ( 1 1 ) )( ( xy xy y x y x r y n x n y x n r n i i n i i n i i i 相関係数 r r r y x y x y x との発現パターンがほぼ正反対との発現パターンがばらばらとの発現パターンが酷似 2) (0 1 D r 距離 D 2 1) ( D r D r D r 15

階層的クラスタリング 1. 遺伝子間距離を計算例 :4 遺伝子の場合距離 D 1 r (0 D 2) 相関係数 r 相関係数 r 相関係数 r... 1,2 1,3 1,4 0.98 0.

16 階層的クラスタリング 1. 遺伝子間距離を計算例 :4 遺伝子の場合距離 D 1 r (0 D 2) 相関係数 r 相関係数 r 相関係数 r... 1,2 1,3 1, 距離 D 1,2 距離 D 距離 D 1 r 距離 D 2 1,3 1,4 (0 D 1) ( 0.01) ( 0.78)

17 階層的クラスタリング 2. 距離行列を作成距離 D 距離 D 距離 D... 1,2 1,3 1, ( 0.01) ( 0.78) 距離行列イメージ 17

18 階層的クラスタリング 3. 樹形図を作成距離行列距離 D D 3,4 二つのクラスター間の距離?!

19 階層的クラスタリング 3. 樹形図を作成平均連結法の場合 D ( D D D ) / 4 4 1, 3 1, 4 2, 3 D2, ( ) / 単連結法の場合 min( D1, 3, D1, 4, D2, 3, D2, 4) 0.47 完全連結法の場合 max( D1, 3, D1, 4, D2, 3, D2, 4)

20 手順 2: サンプル間クラスタリング Affymetrix GeneChip の場合様々な前処理法を適用して得られた遺伝子発現行列データごとに行う結果を眺めて反復実験結果が同一クラスターに含まれる前処理法のデータを採用 RMA がよかった場合 : それだけを採用でよしそれ以外の場合で残りの二つの前処理法の結果のトポロジーが同じ場合 : 二つのデータを同時並行で解析 ( 論文では一つのみ ) 三つの結果がいずれも異なっていた場合 : ご愁傷さまです... その他のメーカーの場合メーカー推奨のやり方に従って遺伝子発現行列データ ( 基本的に一つのみ ) を得る 20

21 Nakai et al., BBB, 2008 クラスタリング結果の解釈例肝臓 (LIV) 白色脂肪 (WAT) 褐色脂肪 (BAT) 通常 (fed) vs. 24 時間絶食 (fasted) MAS-preprocessed data ( メーカー推奨?!) RMA-preprocessed data RMA がいいと判断 ( この場合 ) 21

22 Nakai et al., BBB, 2008 階層的クラスタリング例肝臓 (LIV) 白色脂肪 (WAT) 褐色脂肪 (BAT) 最適なクラスター数 K は? K=2 K=3 K=5 K= K= K=4 K=

23 Ben-Hur et al., PSB, 2002 最適なクラスター数を見積もる方法様々な K について ( 例えば K=2) 全サンプル (n) のクラスタリング結果を K 個に分割した結果とサブサンプル ( 例えば n*0.7) のクラスタリング結果を K 個に分割した結果の類似度を計算全サンプルの結果 1 回目 100 回の結果全て LIV とそれ以外を分割できた場合回目サブサンプリングデータでクラスタリングを例えば 100 回繰り返し 100 回目 23

24 Ben-Hur et al., PSB, 2002 手順 2 : クラスター数をチェック K の値をいくつか試して ( 例では 2~9) 最適な K の値を同定この場合は K=2, 3 が最適なクラスター数言いたいことと同じだったらラッキー 24

25 二群間比較例 1) A 群 : 癌サンプル B 群 : 正常サンプル癌と正常で発現の異なる遺伝子 A A B B x 1, 1 x 1, 2 x 1, 2 x 1, 2 A A B B x 2, 1 x 2, 2 x 2, 2 x 2, 2 A A B B xi, 1 xi, 2 xi, 2 xi, 2 A A B B xn, 1 xn, 2 xn, 2 xn, 2 25

26 Golub et al., Science, 1999 二群間比較解析例 ) 急性白血病 A 群 : リンパ性 (27 サンプル ) B 群 : 骨髄性 (11サンプル) 白血病のタイプで発現の異なる遺伝子群を同定 2009/10/21 26

27 二群間比較 ( 解析手法 ) 倍率変化 (Fold change; FC) に基づくランキング法 2-fold, 3-fold (FC) The limit fold change model (Mutch et al., BMC Bioinformatics, 2002) Rank product (RP; Breitling et al., FEBS Lett., 2004) WAD (Kadota et al., Algorithm. Mol. Biol., 2008) t- 統計量に基づくランキング法 a signal-to-noise statistic (Golub et al., Science, 1999) Student s (or Welch) t-test SAM (samt; Tusher et al., PNAS, 2001) Samroc (Broberg, P., Genome Biol., 2003) a moderated t statistic (Smyth, GK., Stat. Appl. Genet. Mol. Biol., 2004) Intensity-based moderated t statistic (IBMT; Sartor et al., BMC Bioinformatics, 2006) Shrinkage t statistic (Opgen-Rhein and Strimmer, Stat. Appl. Genet. Mol. Biol., 2007) その他 Probability of Positive LogRatio (PPLR; Liu et al., Bioinformatics, 2006) FCPC (Qin et al., Bioinformatics, 2008) 個々の遺伝子の発現変動の度合いを調べる研究 27

28 参考資料二群間比較 (t- 統計量に基づくランキング法 ) 二群間の平均の差が大きく群内のばらつきが小さい遺伝子 iを抽出 a signal-to-noise(s2n) 統計量 R( i) i A U A i B U 二群間の平均の差 A 群内のばらつき B 群内のばらつき対数変換 (log2 変換 ) 後のデータ i B i 標本平均標本分散不偏分散 A j R(1) R(2) R(3) S U A 2 i A i 2 i A 1 n A A n A j1 1 n 1 A n i j 1 na i ( Aj n A 6, n 5, n n 統計量の絶対値が大きい候補発現変動遺伝子 n B A A j1 n A ( A B i j i ) 2 i A ) 2 28

29 参考資料二群間比較 ( 倍率変化に基づくランキング法 ) log 比 :( 対数変換後のデータなので )t 検定系の数式の分子のみに相当 R i i ( i) log( FC) A B 二群間の平均の差対数変換 (log2 変換 ) 後のデータ Average Difference(AD) 統計量と私は呼んでいる R(1) R(2) R(3) 統計量の絶対値が大きい候補発現変動遺伝子 29

Average Difference (AD) 統計量 AD 平均シグナル強度 i Bi A x / 2 i i Bi Ai xを (0~1) の範囲に規格化 xi min( x) wi max( x) min( x) WAD 統計量 WAD i

30 Kadota K, Nakai Y, Shimizu K, AMB., 3:8, 2008 二群間比較 ( 倍率変化に基づくランキング法 ) WAD:log 比を基本としつつ全体的にシグナル強度の高い遺伝子が上位にくるように重みをかけた統計量 unlogged data log 2 -transformed data Average Difference (AD) 統計量 AD 平均シグナル強度 i Bi A x / 2 i i Bi Ai xを (0~1) の範囲に規格化 xi min( x) wi max( x) min( x) WAD 統計量 WAD i 参考資料 AD w i i AD AD i gene6 B i A i (6 7) / 2 ( より 2) / 3 x x i gene6 B A i i (6 7) / 2 (1 2 2) / 3 / / 2より WAD の一位 :gene4, AD の一位 :gene6 w w i xi min( x) より max( x) min( x) gene

31 二群間比較 ( 倍率変化に基づくランキング法 ) Rank products (RP):A 群 vs. B 群の総当たりの比を計算しその順位の相乗平均を統計量とする入力データ総当りの発現比を計算 Breitling et al., FEBS Lett., 2004 (n A n B ) = 9 通り参考資料 n A = 3 n B = 3 列ごとに Rank を計算した後各行に対して相乗平均値 (RPs) を計算 31

32 評価の実際例 :Affymetrix の二群間比較 ( 最もよく研究されている ) Gene Ontology 解析 ( 未知サンプルの ) 分類モチーフ解析パスウェイ解析感度特異度既知の発現変動遺伝子をどれだけ上位にランキング可能か? 再現性同じサンプルの比較結果 ( 発現変動遺伝子リスト ) が場所間でどれだけ一致しているか? 32

33 感度特異度を AUC 値で評価どの前処理法がいい?( 比較例 :MAS5 vs. RMA) 既知の発現変動遺伝子をどれだけ上位にランキング可能か?(AUC 値の高さ ) MAS5 の遺伝子発現行列 log 比を計算 log 比でランキング AUC 値 =100% RMA の遺伝子発現行列 AUC 値 =83.3% 33

34 MAQC Consortium, Nat. Biotechnol., 24: , 2006 再現性を一致度で評価 MicroArray Quality Control (MAQC) プロジェクトで提唱 (0 POG 100%) POG 値が高いランキング結果の頑健性 ( 再現性 ) が高い方法 MAS5 WAD MAS5 WAD 34

35 MAQC Consortium, Nat. Biotechnol., 24: , 2006 POG between 九大 and 東大再現性を一致度で評価 MicroArray Quality Control (MAQC) プロジェクトで提唱 (0 POG 100%) POG 値が高いランキング結果の頑健性 ( 再現性 ) が高い方法上位 x 個の集合 x = 前処理法 :MAS5, ランキング法 :WAD 九大東大前処理法 :MAS5, ランキング法 :samt 九大東大 x 再現性 :WAD > samt 35

36 Kadota K, Nakai Y, Shimizu K, AMB, 4: 7, 2009 再現性解析結果 ( 前処理法 :FARMS) サンプル C 5 例 vs. サンプル D 5 例 Site1 Site4 Site2 上位 100 個の集合 Site5 Site3 Site4 Site6 17% Site1 Site2 x Site5 Site6 Site3 再現性 :WAD > MAQC 推奨法 (AD) 36

37 Kadota K, Nakai Y, Shimizu K, AMB, 4: 7, 2009 結論 (Affymetrix データ ; 二群間比較 ) 感度特異度が高い方法 ( 組合せが重要である!) ( 発現変動遺伝子リストの ) 再現性が高い方法 Fold Change に基づく方法従来 :t- 統計量に基づく方法 ( 前処理法によらず )WAD 従来 : Average Difference (AD) 法 MAQC Consortium, Nat. Biotechnol., 24: , 2006 No Kadota s guidelines, no good research! 37

38 手順 3: 発現変動遺伝子のランキング Affymetrix GeneChip の場合推奨の組み合わせのものを利用 RMA データの場合は Rank products を利用などその他のメーカーの場合 ( 今のところ根拠なし ) チップごとに正規化したデータ WAD 全サンプルのデータを Quantile 正規化したようなデータ Rank products 38

39 クラスタリング結果を眺めることで... 本物 ( 真の発現変動遺伝子 ) があるかどうかの検討がつきます MAS-preprocessed data RMA-preprocessed data 発現変動遺伝子なさそう... 発現変動遺伝子沢山ありそう

40 RMA-quantified data なので... Rank products 法を適用 WAT サンプルの 4 fed vs. 4 fasted samples のデータの解析結果 MAS-preprocessed data RMA-preprocessed data FDR の閾値 0.01 以下 :4 個 (fasted < fed), 45 個 (fasted > fed) 0.10 以下 :90 個 (fasted < fed), 198 個 (fasted > fed) FDR の閾値 0.01 以下 :359 個 (fasted < fed), 278 個 (fasted > fed) 0.10 以下 :970 個 (fasted < fed), 928 個 (fasted > fed)

41 2009/10/21 二群間比較解析戦略発現変動遺伝子 ( マーカー遺伝子 ) の同定個々の遺伝子について統計量を算出しランキング手法選択のガイドライン (Kadota et al., AMB, 2009) 感度特異度重視の場合再現性重視の場合 Gene Set Enrichment Analysis (GSEA) アノテーション情報が豊富な生物種用の解析手段同じセットに属する遺伝子をひとまとめにして解析例 1: 酸化的リン酸化に関係する遺伝子セット (KEGG: hsa00190) 例 2: 脂肪酸 β 酸化に関係する遺伝子セット ( GO: ) 比較する二群間でその遺伝子セットが動いたかどうかを評価帰無仮説 : 動いてない対立仮説 : 動いた沢山の遺伝子セットについて解析を行い動いた遺伝子セットを列挙 positional gene sets pathway gene sets motif gene sets 様々な視点での解析が可能 GO gene sets etc... 41

42 様々な遺伝子セットは MSigDB からゲット例 :KEGG Pathway 遺伝子セット Pathway ID Name Gene symbols 2009/10/21 1 行につき 1 セット 42

43 様々な GSEA 系の解析手法 GSEA (Subramanian et al., PNAS, 2005) PAGE (Kim and Volsky, BMC Bioinformatics, 2005) Hotelling s T 2 -test (Kong et al., Bioinformatics, 2006) GSA (Efron and Tibshirani, Ann. Appl. Stat., 2007) GeneTrail (Backes et al., NAR, 2007) SAM-GS (Dinu et al., BMC Bioinformatics, 2007) GSEA-P (Subramanian et al., Bioinformatics, 2007) GlobalANCOVA (Hummell et al., Bioinformatics, 2008) 2009/10/21 43

44 Kim and Volsky, BMC Bioinformatics, 2005 PAGE 法 Parametric Analysis of Gene set Enrichment の略 1. 各遺伝子 iについて対数変換後のデータのaverage Difference (AD i i i i ) を計算 AD A B ( i 1,2,..., 2. AD i の平均 μ と標準偏差 σ を計算 3. 興味ある遺伝子セット ( 例 :i=5,89, 684, 2543, に相当する計 m 個の遺伝子 ) のADの平均 S m を計算 S m ( AD AD AD AD...) / m 4. Z スコアを計算 Z ( Sm ) m / Z スコアの絶対値が大きい遺伝子セットほど二群間でより発現変動していると解釈, a) 2009/10/21 44

45 a genes GSEA 以前の解析手段例 : 酸化的リン酸化関連遺伝子セット 1. Average Differenceのような統計量を各遺伝子について算出 2. 上位 x 個を抽出し酸化的リン酸化関連遺伝子群のバックグラウンド (b/a) に対する濃縮度合い (c/x) を評価 A 群 B 群 A 群 B 群酸化的リン酸化関連遺伝子 ( チップ中に b 個 ) の位置帰無仮説 : チップ中の全遺伝子数 (a) に対する酸化的リン酸化関連遺伝子数 (b) の割合 (b/a) と酸化的リン酸化関連遺伝子数 (b) に対する上位 x 個の中に占める酸化的リン酸化関連遺伝子数 (c) の割合 (c/x) は等しい 2009/10/21 45

46 a genes GSEA 以前の解析手段の問題点 1 上位 x 個の x 次第で結果が変わる A 群 B 群 A 群 B 群酸化的リン酸化関連遺伝子 ( チップ中に b 個 ) の位置 2009/10/21 46

47 a genes GSEA 以前の解析手段の問題点 2 下図のように全体としては酸化的リン酸化関連遺伝子セットが有意差があるといえるような場合でも上位 x 個の中に一つも含まれないので有意差があるといえなくなる現実の解析では酸化的リン酸化関連遺伝子セットが動いていることを見落とす A 群 B 群 A 群 B 群酸化的リン酸化関連遺伝子 ( チップ中に b 個 ) の位置 2009/10/21 47

48 様々な GSEA 系手法なぜ次々と提案されるのか? Ans.1: 発現変動遺伝子のランキング法 (gene-level statistics) はいくらでもある PAGE:Average Difference (AD) 倍率変化そのもの GSEA:S2N 統計量などその他 :Rank products, WAD, SAM など Ans.2: 興味ある遺伝子セットの偏り度合い ( 濃縮度 ) を見積もる統計量 (gene set statistics) はいくらでもある PAGE:Z 検定 GSEA:Enrichment Score その他 : 平均 % 順位, AUC, median など Ans.3: 有意性を評価する手段もいくつか考えられる sample label permutation gene resampling 極論 : 論文になっていない組合せを新規手法だ! とすることも可能 /10/21 48

49 Ackermann and Strimmer, BMC Bioinformatics, 2009 手法選択のガイドラインはない ( に等しい ) どの遺伝子セットが動いているいないという正解情報 ( 地上の真実 ) を知るすべがない論文でありがちなプレゼンテーション既知の遺伝子セットはちゃんと上位にあった我々はさらに他に動いている遺伝子セットを見つけた ( 感度の高さをアピール ) 感度の高さという点については正しいのかもしれないが特異度は低いのかも... ( 本当は動いていない遺伝子セットまで動いていると判断してしまうこと ) シミュレーションで本当は動いていないデータセットを作成することはできるがその結果と現実の結果には相当のギャップがある 2009/10/21 49

50 GSEA 系手法を使えるのはごく一部の生物種アノテーション情報が豊富な生物種は Gene Ontology やパスウェイの情報が豊富多くの遺伝子セットを用意できる GSEA 系手法を適用可能それ以外の生物種はまずは様々な発現変動遺伝子をひたすら同定しまくるなどして地道にアノテーション情報を増やしていく以外にない ( のではないだろうか ) 2009/10/21 50

51 手順 4-1:GSEA を実行重複した Gene symbol 名のものをまとめたファイルを作成 31,099 行 (data_rma.txt) 14,140 行 (data_rma_nr.txt) 理由 1: 変なバイアスを除きたいから理由 2: 遺伝子セットが Gene symbol で与えられているから 51

52 手順 4-1:GSEA を実行必要なファイルを MSigDB からダウンロード 52

53 GSEA 実行例 ( 手順 4-2) LIV サンプルの 4 fed vs. 4 fasted samples データの Gene Ontology(Biological Process) 解析結果上位 10 遺伝子セット絶対値の大きいものほど偏り度合いが高いことを表す符号はその方向正の値 A 群 > B 群 p 値論文の表の完成この GO ID に含まれる遺伝子セットのメンバー数 63 個中 52 個が自分が用いたアレイ中に搭載されている 53

54 手順 4-3:GO の階層構造にマップ LIV サンプルの 4 fed vs. 4 fasted samples データの Gene Ontology (Biological Process) 解析結果上位 x 遺伝子セット ( 例 :x = 10) の GO IDs を QuickGO にかける 54

55 手順 4-3:GO の階層構造にマップ LIV サンプルの 4 fed vs. 4 fasted samples データの Gene Ontology (Biological Process) 解析結果上位 x 遺伝子セット ( 例 :x = 10) の GO IDs を QuickGO にかける論文の図の完成 55

56 GSEA 実行例 ( 手順 4-2 ) LIV サンプルの 4 fed vs. 4 fasted samples データの KEGG Pathway 解析結果上位 10 遺伝子セット論文の表の完成 56

57 手順 4-3 : パスウェイ上にマップ LIV サンプルの 4 fed vs. 4 fasted samples データの KEGG Pathway 解析結果から第 3 位の HSA00071 を構成する遺伝子メンバーの二群間 (fed vs. fasted) での変動の程度を 4 諧調色で表示 logratio <= -1 を水色 (A 群で発現上昇 ) -1 < logratio < 0 を薄水色 (A 群で発現上昇 ) 0 < logratio < 1 を薄ピンク色 (B 群で発現上昇 ) logratio >= 1 をピンク色 (B 群で発現上昇 ) logratio = mean(b) mean(a) 57

58 水色 :A 群で発現上昇桃色 :B 群で発現上昇 58

59 問題点 EC 番号と Gene symbol が 1 対 1 対応ではない... 例 ) HSA00071 を構成する 39 gene symbols のうち EC: に対応するのは 4 つある... 現状では最終的に反映されている色は同一 EC 番号の一番最後に出てきた gene symbol (i.e., CPT2) の発現レベル 59

60 アグリバイオインフォマティクス教育研究プログラムのフォーラム活動について本プログラムでは研究課題ごとにフォーラムを形成しセミナーシンポジウムの開催から企業との共同研究学位論文の指導などを行い当該課題の研究教育の活性化を図りますフォーラムのメンバーは本研究科の教員のほか他大学企業試験研究機関の方々から構成されますこれらのメンバーから農学生命情報科学実習 II の受講を通して学位論文の研究におけるバイオインフォマティクスに関係した研究の指導を受けることができますバイオインフォマティクスを利用した農学生命科学の研究あるいはバイオインフォマティクスそのものの研究を行って学位を取得した人には修了認定証を発行します修了の認定は各専攻の学位審査とは別にフォーラムのメンバーが審査会を開いて行います研究指導は研究室の指導教員との合意に基づいて行いますので希望する人は指導教員と相談の上アグリバイオインフォマティクス教育研究プログラム事務局までご連絡下さい現在のところ以下の 4 つのフォーラムが形成されています : 微生物インフォマティクスフォーラム基盤バイオインフォマティクスフォーラムアグリ / バイオセンシングと空間情報学フォーラム食品インフォマティクスフォーラム 2009/10/21 60

61 遺伝子発現行列二群間比較様々な組織 ( 条件 ) 時系列データ A A B B x 1, 1 x 1, 2 x 1, 2 x 1, 2 A A B B x 2, 1 x 2, 2 x 2, 2 x 2, 2 x 1,1 x 2,1 x 1,2 x 2,2 x 1,3 x 2,3 x 1,4 x 2,4 x 1,1 x 2,1 x 1,2 x 2,2 x 1,3 x 2,3 x 1,4 x 2,4 A xi, 1 A B B xi, 2 xi, 2 xi, 2 x i,1 x i,2 x i,3 x i,4 x i,1 x i,2 x i,3 x i,4 A A B B xn, 1 xn, 2 xn, 2 xn, 2 x n,1 x n,2 x n,3 x n,4 x n,1 x n,2 x n,3 x n,4 61

62 組織特異的遺伝子検出法ランキングに基づく方法 Dixon test (Greller and Tobin, Genome Res., 9, , 1999) Pattern matching(pavlidis and Noble, Genome Biol., 2, research0042, 2001) Entropy (Schug et al., Genome Biol., 6, R33, 2005) Tissue specificity Index (Yanai et al., Bioinformatics, 21, , 2005) 外れ値検出に基づく方法 Akaike s Information Criterion (AIC) (Kadota et al., Physiol. Genomics, 12, , 2003) Sprent s non-parametric method (Ge et al., Genomics, 86, , 2005) その他 Tukey-Kramer s Honest Significance Difference (HSD) test (Liang et al., Physiol. Genomics, 26, , 2006) ROKU (Kadota et al., BMC Bioinformatics, 7, 294, 2006) 62

63 組織特異的遺伝子検出法様々な前処理 ( 正規化 ) 法様々な二群間での発現変動遺伝子検出法重視すべき評価基準は? 感度特異度再現性推奨ガイドライン結論 : おすすめは ROKU 63

64 心臓胃大脳心臓胃大脳心臓胃大脳組織特異的遺伝子検出法やりたいこと 1 x 1,1 x 1,2 x 1,3 x 1,4 大脳特異的高発現遺伝子 x 2,1 x 2,2 x 2,3 x 2,4 x i,1 x i,2 x i,3 x i,4 x n,1 x n,2 x n,3 x n,4 心臓と大脳特異的高発現遺伝子入力 : 遺伝子発現行列出力 : 任意の組織特異的遺伝子様々な特異的発現パターンを組織特異性の度合いで統一的にランキングしたい 64

al., Genome Biol., 2005 H( x) 0 H( x) 1.

65 組織特異的遺伝子検出法 2 エントロピーによるランキング遺伝子 x = (x 1, x 2,, x n ) のエントロピー H(x) H x ) p log ( p ), where p x H(x) のとりうる範囲 : 0 H(x) log 2 (n) n ( i i i i i x 1 2 i Schug et al., Genome Biol., 2005 H( x) 0 H( x) H( x) 1.45 H ( x) H ( x ) 3.32 log 2( n ) エントロピーが低い組織特異性が高いエントロピーが高い組織特異性が低いエントロピーでランキングすることにより複数外れ値に対応可能 65

66 Schug et al., Genome Biol., エントロピー計算例遺伝子 i のエントロピー H(x i ) H( x i N ) 1 pij log 2( p 0 H log 2 N j ij ) p ij x ij / N j 1 x ij 0 H 2.32 特異的発現パターン低いエントロピーそうでないパターン高いエントロピー 66

67 Schug et al., Genome Biol., 2005 参考資料組織特異的遺伝子検出法 2 エントロピーの短所 1. 組織特異的低発現パターンなどの検出が不可能 0 H(x) log 2 (n) 3.32 H( x) 3.29 H( x) H( x) 特異的組織の同定が不可能上位にランキングされない H( x) 0 H( x) 0 H( x) 0 どの組織で特異的なのか分からない 67

68 Kadota et al., BMC Bioinformatics, 2006 参考資料組織特異的遺伝子検出法 3ROKU 1. 遺伝子発現ベクトル x を変換 : x x by x i = x i T bw 0 H(x) log 2 (n) 3.32 H( x ) 1.48 H( x ) H( x ) AIC に基づく外れ値検出法を採用上位にランキングされるどの組織で特異的なのか分かる 68

AIC nn log 2 no nn 入力 n n n n o n n o : Outlier ( 外れ値 ) の数 : Non - outlier の数

69 組織特異的遺伝子検出法 4AIC に基づく外れ値検出法 Akaike s Information Criterion (AIC) 様々な外れ値の組み合わせモデルからAICが最小の組み合わせ (MAICE) を探索計算例 : log nn! AIC nn log 2 no nn 入力 n n n n o n n o : Outlier ( 外れ値 ) の数 : Non - outlier の数 ˆ : 標準偏差 ( n) : サンプル数出力上田太一郎, 応用統計学, 1996 Kadota et al., Physiol. Genomics, 2003 低発現側の外れ値 :-1, 高発現の ~:1, それ以外 :0 69

70 参考資料組織特異的遺伝子検出法 4AIC に基づく外れ値検出法様々な外れ値の組み合わせモデルから AIC が最小の組み合わせ (MAICE) を探索様々な外れ値の組み合わせモデル最大探索範囲 Nmax = n/2 = 5 AIC n n n n o n n n n o log ˆ : 標準偏差 2 n ( n) : サンプル数 : Outlier ( 外れ値 ) の数 : Non - outlier の数 o log n n n n! 70

機能ゲノム学（第6回）

機能ゲノム学（第6回）トランスクリプトーム解析手法の開発東京大学大学院農学生命科学研究科アグリバイオインフォマティクス人材養成ユニット門田幸二 2008/12/08 トランスクリプトーム (transcrptome) とは細胞中に存在する転写物全体 (transcrpt + ome) トランスクリプトーム解析技術 DNA マイクロアレイ Affymetrx GeneChp, cdna アレイ, 電気泳動に基づく方法