農学生命情報科学特論I

Transcription

1 版ほぼ完成です 最終回のレポート課題はありません 基本的な考え方と解析戦略の変遷 ( スライド 3-29 あたり ) は確実に省略しますので講義前に 3-29 については自分で見ておいてください スライド についても残り時間次第です 最終回ですので アンケートのほうもよろしくお願いします 農学生命情報科学特論 I 第 4 回 大学院農学生命科学研究科アグリバイオインフォマティクス教育研究プログラム 2 微生物科学イノベーション連携研究機構門田幸二 ( かどたこうじ ) kadota@iu.a.u-tokyo.ac.jp July 3, 208

2 講義予定 第 回 (208 年 06 月 2 日 ) カウント情報取得の続き データの正規化 (RPK, RPM, RPKM/FPKM) 第 2 回 (208 年 06 月 9 日 ) サンプル間クラスタリング R のクラスオブジェクト R の Reference Manual の読み解き方 クラスタリング結果の客観的な評価 第 3 回 (208 年 06 月 26 日 ) 発現変動解析 ( 反復あり / なしの 2 群間比較 ) M-A plot 発現変動解析 :3 群間比較 デザイン行列 第 4 回 (208 年 07 月 03 日 ) 機能解析 ( 発現変動遺伝子セット解析 ) GSEA MSigDB GSVA の実行 July 3, 208 2

3 Contents 機能解析 ( 発現変動遺伝子セット解析 ) 全体像 基本的な考え方と解析戦略の変遷 様々なプログラム 遺伝子セット情報の取得 (gmtファイルの取得) 発現データ情報と遺伝子セット情報のIDの対応付け 検証用 RNA-seqカウントデータセットPickrell data( なぜGSVAにしたか ) GSVAの解説 PDFを読み解く ( 手元のc.all.v6..entrez.gmt をどう読み込ませるか ) GSVAdataパッケージ提供の MSigDB c2コレクションであるc2broadsetsを理解する 手元のgmtファイルを読み込ませて GeneSetCollection 形式で取り扱えるようにする GSVA の解説 PDF を読み解く ( 手元の発現データファイルをどう取り扱うか ) ExpressionSetの取り扱い nsfilter 関数を用いた同一 IDの重複除去 メインプログラムgsva 関数が入力として受け付けるデータ形式 (ExpressionSetとMatrix) 検証用 RNA-seqカウントデータセットPickrell dataのイントロ スルーしていいところ MSigDB c2コレクションに2つの性特異的遺伝子セットを追加したものでgsvaを実行 ユニークなEntrez gene IDで グループごとに分離させた発現データファイル作成 整形後の発現データファイルとc.all.v6..entrez.gmtを入力としてGSVAを実行 July 3, 208 3

4 機能解析 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 をもつ遺伝子セットが G 群 vs. G2 群比較で動いていた July 3, 208 4

5 機能解析の全体像 で機能解析 ( 発現変動遺伝子セット解析 ) に関する全体像を述べています July 3, 208 5

6 機能解析の全体像 第 世代の機能解析 (ORA) 2GSEA が含まれる第 2 世代の機能解析 (FCS) ORA や FCS の用語を含む詳細については後述 と書いてありますが講義では省きます 2 July 3, 208 6

7 遺伝子セット次第で ここがキモ 2 遺伝子セット次第で GO 解析や Pathway 解析にもなり得る 2 July 3, 208 7

8 MSigDB MSigDB というサイトで 遺伝子セット情報の.gmt という拡張子のついたファイルが提供されています July 3, 208 8

9 MSigDB MSigDB というサイトで 遺伝子セット情報の.gmt という拡張子のついたファイルが提供されています July 3, 208 9

10 gmt ファイルを入手 MSigDB というサイトで 遺伝子セット情報の.gmt という拡張子のついたファイルを予め入手しておかねばならない July 3, 208 0

11 入力ファイルは 2 つ 発現変動遺伝子検出の場合は 入力が つ ( 発現行列データ ) であった 発現変動遺伝子セット解析の場合は 発現データファイルに加えて どの遺伝子がどの遺伝子セットに属するかという情報を含む gmt ファイルも必要です July 3, 208

12 様々なプログラムがある 2005 年の GSEA 論文発表前後を含めて 様々な方法が提案されています July 3, 208 2

13 Contents 機能解析 ( 発現変動遺伝子セット解析 ) 全体像 基本的な考え方と解析戦略の変遷 様々なプログラム 遺伝子セット情報の取得 (gmtファイルの取得) 発現データ情報と遺伝子セット情報のIDの対応付け 検証用 RNA-seqカウントデータセットPickrell data( なぜGSVAにしたか ) GSVAの解説 PDFを読み解く ( 手元のc.all.v6..entrez.gmt をどう読み込ませるか ) GSVAdataパッケージ提供の MSigDB c2コレクションであるc2broadsetsを理解する 手元のgmtファイルを読み込ませて GeneSetCollection 形式で取り扱えるようにする GSVA の解説 PDF を読み解く ( 手元の発現データファイルをどう取り扱うか ) ExpressionSetの取り扱い nsfilter 関数を用いた同一 IDの重複除去 メインプログラムgsva 関数が入力として受け付けるデータ形式 (ExpressionSetとMatrix) 検証用 RNA-seqカウントデータセットPickrell dataのイントロ スルーしていいところ MSigDB c2コレクションに2つの性特異的遺伝子セットを追加したものでgsvaを実行 ユニークなEntrez gene IDで グループごとに分離させた発現データファイル作成 整形後の発現データファイルとc.all.v6..entrez.gmtを入力としてGSVAを実行 July 3, 208 3

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

15 N =0,000 genes 基本的な考え方 遺伝子ごとの発現変動の度合いを数値化 例 :t- 統計量 log2(g2/g) 相関係数 基本はデフォルトでやるようですが 様々な選択肢があります G 群 G2 群 G 群 G2 群 発現変動遺伝子 (G >> G2) 変動していない遺伝子 (G == G2) 発現変動遺伝子 (G << G2) July 3, 208 5

16 N =0,000 genes 基本的な考え方 発現変動順にソート後の酸化的リン酸化関連遺伝子セットのステレオタイプな分布 どうやって偏りを評価するのか? 変動している 変動していない G 群 G2 群 G 群 G2 群 July 3, 208 6

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

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

19 rcode_ora_basic.txt ORA( 超幾何検定 ) DEG として 500 個抽出したとき 酸化的リン酸化関連遺伝子が 6 個以上含まれる確率として算出 N=0000 個の遺伝子発現データ中に XXX= 酸化的リン酸化関連遺伝子は n=7 個含まれていた 上位 X=500 個の発現変動遺伝子 (DEG) の中に x=6 個の酸化的リン酸化関連遺伝子が含まれていた 帰無仮説 : 酸化的リン酸化関連遺伝子の割合は DEG と non-deg 間で差がない ORA として Fisher s hypergeometric test を利用 (Tavazoie et al., Nat Genet., 22: , 999) July 3, 208 9

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

21 rcode_ora_basic.txt ORA( カイ二乗検定 ) DEG として 500 個抽出したとき 酸化的リン酸化関連遺伝子が 6 個以上含まれる確率として算出 July 3, 208 2

22 rcode_ora_basic.txt 直感は重要 N=0000 個の遺伝子発現データ中に XXX= 酸化的リン酸化関連遺伝子は n=7 個存在する 2 上位 X=500 個の発現変動遺伝子 (DEG) の中に x=6 個の酸化的リン酸化関連遺伝子が含まれていた 3 帰無仮説 : 酸化的リン酸化関連遺伝子の割合は DEG と non-deg 間で差がない の段階で 調べたい遺伝子セットは 7/0,000 = 0.07% の割合だと考える 2 で 6/,500 = 0.4% の割合に濃縮されていると考える 3 今やっているのは ある 2 群間比較 もし比較している群間で この遺伝子セットが全体として発現変動していなかったとしたら ランダムで,500 個とった時に この遺伝子セット中の遺伝子が含まれる割合は 0.07% なので 個数だと % =.05 個程度しか含まれないはず 実際に得られたのは 6 個なので 偶然こんな結果が得られたとは考えにくい 起こるとしたら 4 くらい低い確率なんだね だから 発現変動遺伝子セット と考えよう という思考回路 4 July 3,

23 N =0,000 genes ORA Over-Representation Analysis (ORA) 上位 500 個のうち 酸化的リン酸化関連遺伝子が 7 個中 4 つ以上含まれていれば p < 0.05 で検出可能ということを意味する G 群 G2 群 G 群 G2 群 p < 0.05 を灰色で示した July 3,

24 ORA Over-Representation Analysis (ORA) GenMAPP (Dahlquist et al., Nature Genet., 3: 9-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, 200) GenMAPP は比較的有名 July 3,

25 N =0,000 genes 第 世代 (ORA) の短所 全体的には動いているものの 個々の発現変動の度合いが弱い場合に検出困難 2 上位 X 個の X 次第で結果が変わる 3 情報量低下 ( 発現変動の度合い カウント情報 ) 3 G 群 G2 群 2 G 群 G2 群 July 3,

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

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

28 第 2 世代 (FCS) Functional Class Scoring (FCS) GSEA (Subramanian et al., PNAS, 02: , 2005) PAGE (Kim and Volsky, BMC Bioinformatics, 6: 44, 2005) sigpathway (Tian et al., PNAS, 02: , 2005) GSA (Efron and Tibshirani, Ann. Appl. Stat., : 07-29, 2007) GeneTrail (Backes et al., Nucleic Acids Res., 35: W86-W92, 2007) SAM-GS (Dinu et al., BMC Bioinformatics, 8: 242, 2007) 最も有名なのが GSEA ここでリストアップされているのは 基本的にマイクロアレイデータ解析用なので情報としては古い よって (R で ) 塩基配列解析ではほとんどリストアップしていない July 3,

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

30 Contents 機能解析 ( 発現変動遺伝子セット解析 ) 全体像 基本的な考え方と解析戦略の変遷 様々なプログラム 遺伝子セット情報の取得 (gmtファイルの取得) 発現データ情報と遺伝子セット情報のIDの対応付け 検証用 RNA-seqカウントデータセットPickrell data( なぜGSVAにしたか ) GSVAの解説 PDFを読み解く ( 手元のc.all.v6..entrez.gmt をどう読み込ませるか ) GSVAdataパッケージ提供の MSigDB c2コレクションであるc2broadsetsを理解する 手元のgmtファイルを読み込ませて GeneSetCollection 形式で取り扱えるようにする GSVA の解説 PDF を読み解く ( 手元の発現データファイルをどう取り扱うか ) ExpressionSetの取り扱い nsfilter 関数を用いた同一 IDの重複除去 メインプログラムgsva 関数が入力として受け付けるデータ形式 (ExpressionSetとMatrix) 検証用 RNA-seqカウントデータセットPickrell dataのイントロ スルーしていいところ MSigDB c2コレクションに2つの性特異的遺伝子セットを追加したものでgsvaを実行 ユニークなEntrez gene IDで グループごとに分離させた発現データファイル作成 整形後の発現データファイルとc.all.v6..entrez.gmtを入力としてGSVAを実行 July 3,

31 GO 解析用 遺伝子セットとして GO の情報を用いる場合 July 3, 208 3

32 GO 解析用 SeqGSA はマニュアルがないに等しい 2GSAR は 発現変動遺伝子セットを探すというよりは 興味ある遺伝子セットを与えてネットワーク図を描き どの遺伝子が hub(hub genes) かを返すのがメイン の講義資料でも解説しています でページ内検索してください 2 July 3,

33 EGSEA EGSEA は 2 複数のツールを組み合わせるやり方 様々な分野でこの種の戦略がよいことは実証されており おそらく妥当 しかしその分だけ依存関係が複雑になるため 私もまだ試してはいない 2 July 3,

34 SeqGSEA SeqGSEA は 手順が煩雑な上 ものすごく計算時間がかかる の講習会資料作成当時の個人の感想です July 3,

35 GSVA GSVA は 2EGSEA でも利用されている また 3 引用回数も多い ( 使いやすい ) ので 後半はこれをベースに説明 2 July 3,

36 パスウェイ解析用 遺伝子セットとして KEGG Pathway や Reactome などの情報を用いる場合 July 3,

37 パスウェイ解析用 GO 解析もできるのに ここにしか記載していないものもいくつかあるはずですのでご注意ください July 3,

38 Contents 機能解析 ( 発現変動遺伝子セット解析 ) 全体像 基本的な考え方と解析戦略の変遷 様々なプログラム 遺伝子セット情報の取得 (gmtファイルの取得) 発現データ情報と遺伝子セット情報のIDの対応付け 検証用 RNA-seqカウントデータセットPickrell data( なぜGSVAにしたか ) GSVAの解説 PDFを読み解く ( 手元のc.all.v6..entrez.gmt をどう読み込ませるか ) GSVAdataパッケージ提供の MSigDB c2コレクションであるc2broadsetsを理解する 手元のgmtファイルを読み込ませて GeneSetCollection 形式で取り扱えるようにする GSVA の解説 PDF を読み解く ( 手元の発現データファイルをどう取り扱うか ) ExpressionSetの取り扱い nsfilter 関数を用いた同一 IDの重複除去 メインプログラムgsva 関数が入力として受け付けるデータ形式 (ExpressionSetとMatrix) 検証用 RNA-seqカウントデータセットPickrell dataのイントロ スルーしていいところ MSigDB c2コレクションに2つの性特異的遺伝子セットを追加したものでgsvaを実行 ユニークなEntrez gene IDで グループごとに分離させた発現データファイル作成 整形後の発現データファイルとc.all.v6..entrez.gmtを入力としてGSVAを実行 July 3,

39 遺伝子セット情報取得 遺伝子セット情報は GMT という形式 ( 拡張子が.gmt) で提供されており 23 つの手段で取得可能です 2 July 3,

40 MSigDB ver. 6. 最も有名なのは MSigDB 2208 年 6 月現在の ver. 6. では 8 個の主要なコレクションが提供されている 2 July 3,

41 MSigDB ver. 6.. H: hallmark gene sets (50 gene sets) 2. c: positional gene sets (326 gene sets) ヒト染色体の位置ごとの遺伝子セットリストファイル (326 gene sets) 3. c2: curated gene sets (4,738 gene sets) CGP: chemical and genetic perturbations (3,409 gene sets) CP: canonical pathways (,329 gene sets) CP:BIOCARTA: BioCarta gene sets (27 gene sets) CP:KEGG: KEGG gene sets (86 gene sets) CP:REACTOME: Reactome gene sets (674 gene sets) 4. c3: motif gene sets (836 gene sets) MIR: microrna targets (22 gene sets) TFT: transcription factor targets (65 gene sets) 5. c4: computational gene sets (858 gene sets) CGM: cancer gene neighborhoods (427 gene sets) CM: cancer modules (43 gene sets) 6. c5: gene ontology (GO) gene sets (5,97 gene sets) BP: biological process (4,436 gene sets) CC: cellular component (580 gene sets) MF: molecular function (90 gene sets) 7. c6: oncogenic signatures gene sets (89 gene sets) 8. c7: immunologic signatures gene sets (4,872 gene sets) 発現変動と関連する KEGG パスウェイを調べたいときは 3 番目の c2 というカテゴリーに属する CP:KEGG というところの 86 個の遺伝子セットが含まれる.gmt ファイルを予めダウンロードしておく July 3, 208 Subramanian et al., PNAS, 02: ,

42 MSigDB ver. 6.. H: hallmark gene sets (50 gene sets) 2. c: positional gene sets (326 gene sets) ヒト染色体の位置ごとの遺伝子セットリストファイル (326 gene sets) 3. c2: curated gene sets (4,738 gene sets) CGP: chemical and genetic perturbations (3,409 gene sets) CP: canonical pathways (,329 gene sets) CP:BIOCARTA: BioCarta gene sets (27 gene sets) CP:KEGG: KEGG gene sets (86 gene sets) CP:REACTOME: Reactome gene sets (674 gene sets) 4. c3: motif gene sets (836 gene sets) MIR: microrna targets (22 gene sets) TFT: transcription factor targets (65 gene sets) 5. c4: computational gene sets (858 gene sets) CGM: cancer gene neighborhoods (427 gene sets) CM: cancer modules (43 gene sets) 6. c5: gene ontology (GO) gene sets (5,97 gene sets) BP: biological process (4,436 gene sets) CC: cellular component (580 gene sets) MF: molecular function (90 gene sets) 7. c6: oncogenic signatures gene sets (89 gene sets) 8. c7: immunologic signatures gene sets (4,872 gene sets) 発現変動と関連する GO の biological process を調べたいときは 6 番目の c5 というカテゴリーに属する BP というところの 4,436 個の遺伝子セットが含まれる.gmt ファイルを予めダウンロードしておく July 3,

43 MSigDB ver. 6.. H: hallmark gene sets (50 gene sets) 2. c: positional gene sets (326 gene sets) ヒト染色体の位置ごとの遺伝子セットリストファイル (326 gene sets) 3. c2: curated gene sets (4,738 gene sets) CGP: chemical and genetic perturbations (3,409 gene sets) CP: canonical pathways (,329 gene sets) CP:BIOCARTA: BioCarta gene sets (27 gene sets) CP:KEGG: KEGG gene sets (86 gene sets) CP:REACTOME: Reactome gene sets (674 gene sets) 4. c3: motif gene sets (836 gene sets) MIR: microrna targets (22 gene sets) TFT: transcription factor targets (65 gene sets) 5. c4: computational gene sets (858 gene sets) CGM: cancer gene neighborhoods (427 gene sets) CM: cancer modules (43 gene sets) 6. c5: gene ontology (GO) gene sets (5,97 gene sets) BP: biological process (4,436 gene sets) CC: cellular component (580 gene sets) MF: molecular function (90 gene sets) 7. c6: oncogenic signatures gene sets (89 gene sets) 8. c7: immunologic signatures gene sets (4,872 gene sets) 講義で利用する解析プログラム GSVA の検証用としては 326 遺伝子セットが最適なので これをダウンロードしておきます July 3,

44 ダウンロード MSigDB 本家からもダウンロードできますが 一気にやると迷惑をかけるので 2 から c.all.v6..entrez.gmt をダウンロードしておいてください Entrez gene ID のヒト染色体の位置ごとの遺伝子セットリストファイル (326 gene sets) です 3 c.all.v6..symbols.gmt もダウンロードしておきましょう 3 2 July 3,

45 c.all.v6..entrez.gmt c.all.v6..entrez.gmt を Excel で眺めるとこんな感じ GMT 形式は 2 列目が gene set 名 32 列目が情報源の URL そして 43 列目以降が遺伝子セットを構成する Entrez gene IDs July 3,

46 c.all.v6..entrez.gmt chr5q23 が 番目のデータセット 2 赤枠部分が chr5q23 という遺伝子セットに含まれる遺伝子群の Entrez gene IDs 2 July 3,

47 chr5q23 最後の Entrez gene ID である 272 まで ずーっと右のほうに移動したところ 遺伝子セットによって 構成メンバー数が異なることがわかります 22 番目の遺伝子セット (chr6q24) は 番目の遺伝子の遺伝子セット (chr5q23) よりも構成メンバー数が多いことがわかります 2 July 3,

48 chr5q23 2 列目の情報源の URL を眺めると July 3,

49 chr5q23 2 列目の情報源の URL を眺めると 2 こんな感じで詳細情報が得られます 3 遺伝子セット chr5q23 のメンバー数は 84 genes であることがわかります 2 3 July 3,

50 Show members で をクリックすると July 3,

51 Show members で こんな感じになって Entrez gene ID に対応する 2gene symbol や 3gene description が見られます 2 3 July 3, 208 5

52 5759 は PTMAP2 番目のデータセット chr5q23 に含まれるメンバーで Entrez gene ID が の 3gene symbol は PTMAP2 であることがわかります July 3,

53 5759 は PTMAP2 参考 c.all.v6..symbols.gmt だとこんな感じになります このファイルは遺伝子セット情報を gene symbols で提供しているものなので妥当ですね July 3,

54 参考 326 gene sets なので 326 行目までで終わりです July 3,

55 326 gene sets なので c.all.v6..entrez.gmt でも当然同じです July 3,

56 Contents 機能解析 ( 発現変動遺伝子セット解析 ) 全体像 基本的な考え方と解析戦略の変遷 様々なプログラム 遺伝子セット情報の取得 (gmtファイルの取得) 発現データ情報と遺伝子セット情報のIDの対応付け 検証用 RNA-seqカウントデータセットPickrell data( なぜGSVAにしたか ) GSVAの解説 PDFを読み解く ( 手元のc.all.v6..entrez.gmt をどう読み込ませるか ) GSVAdataパッケージ提供の MSigDB c2コレクションであるc2broadsetsを理解する 手元のgmtファイルを読み込ませて GeneSetCollection 形式で取り扱えるようにする GSVA の解説 PDF を読み解く ( 手元の発現データファイルをどう取り扱うか ) ExpressionSetの取り扱い nsfilter 関数を用いた同一 IDの重複除去 メインプログラムgsva 関数が入力として受け付けるデータ形式 (ExpressionSetとMatrix) 検証用 RNA-seqカウントデータセットPickrell dataのイントロ スルーしていいところ MSigDB c2コレクションに2つの性特異的遺伝子セットを追加したものでgsvaを実行 ユニークなEntrez gene IDで グループごとに分離させた発現データファイル作成 整形後の発現データファイルとc.all.v6..entrez.gmtを入力としてGSVAを実行 July 3,

57 MSigDB 提供ファイル MSigDB は Entrez gene IDs(.entrez.gmt) と Gene symbols(.symbols.gmt) の 2 種類を提供 当然 本家もそうなっています July 3,

58 機能解析の全体像 機能解析の全体像に関する説明では 2 のあたりで記載しています 2 July 3,

59 機能解析の全体像 手元の発現データの gene ID の種類に応じて 利用する遺伝子セットファイル (.entrez.gmt or.symbols.gmt ) を切り替えます この後で実行する GSVA パッケージの例題の発現データは Entrez gene ID なので c.all.v6..entrez.gmt をダウンロードしたのです July 3,

60 参考 gene symbols の場合 手元の発現データが Entrez gene ID 以外であり 例えば Ensembl gene ID だった場合は Ensembl gene ID と gene symbols の対応情報を取得しておきます そして 発現データのほうを gene symbols に変換しておいて.symbols.gmt ファイルを用いて遺伝子セット解析を行います July 3,

61 参考 gene symbols の場合 様々な ID の変換を行う際によく用いられる BioMart というデータベースを R 経由で実行可能な biomart という R パッケージもあります が話がそれていくので深入りしません July 3, 208 6

62 Contents 機能解析 ( 発現変動遺伝子セット解析 ) 全体像 基本的な考え方と解析戦略の変遷 様々なプログラム 遺伝子セット情報の取得 (gmtファイルの取得) 発現データ情報と遺伝子セット情報のIDの対応付け 検証用 RNA-seqカウントデータセットPickrell data( なぜGSVAにしたか ) GSVAの解説 PDFを読み解く ( 手元のc.all.v6..entrez.gmt をどう読み込ませるか ) GSVAdataパッケージ提供の MSigDB c2コレクションであるc2broadsetsを理解する 手元のgmtファイルを読み込ませて GeneSetCollection 形式で取り扱えるようにする GSVA の解説 PDF を読み解く ( 手元の発現データファイルをどう取り扱うか ) ExpressionSetの取り扱い nsfilter 関数を用いた同一 IDの重複除去 メインプログラムgsva 関数が入力として受け付けるデータ形式 (ExpressionSetとMatrix) 検証用 RNA-seqカウントデータセットPickrell dataのイントロ スルーしていいところ MSigDB c2コレクションに2つの性特異的遺伝子セットを追加したものでgsvaを実行 ユニークなEntrez gene IDで グループごとに分離させた発現データファイル作成 整形後の発現データファイルとc.all.v6..entrez.gmtを入力としてGSVAを実行 July 3,

63 検証用データ 以前のスライドで遺伝子セット解析を実行するプログラムとして 2GSVA を採択した根拠を示した際に 赤枠内のプログラムについては言及したが 3 については言及しなかった 2 3 July 3,

64 Pickrell data 3 2 講義資料作成当時の私は 説明しやすい検証用データセットも探していた AbsFilterGSEA 論文中の 2 のあたりで 3Pickrell data という 2 群間比較用 (29 males vs. 40 females) のカウントデータが存在することを発見 July 3, 208 AbsFilterGSEA(Yoon et al., PLoS One, : e06599, 206) 64

65 MSigDB C この Pickrell data と 2MSigDB C カテゴリーの発現変動遺伝子セット解析結果として 3chryq がおそらく最上位に近い有意性を示し 4male 群で高発現であることまで分かった この段階で c.all.v6..entrez.gmt または c.all.v6..symbols.gmt が有望株だと認識 July 3, 208 AbsFilterGSEA(Yoon et al., PLoS One, : e06599, 206) 65

66 AbsFilterGSEA AbsFilterGSEA が使いやすければよいが July 3, 208 AbsFilterGSEA(Yoon et al., PLoS One, : e06599, 206) 66

67 AbsFilterGSEA 一般に CRAN 提供パッケージは Bioconductor 提供パッケージに比べて利用法の解読が難解です 実際私は 2 の Reference manual をみてガッカリし 他にわかりやすいパッケージはないか探した結果として GSVA を眺め 採用しました 2 July 3, 208 AbsFilterGSEA(Yoon et al., PLoS One, : e06599, 206) 67

68 GSVA 2GSVA は引用回数も多く 3Bioconductor のパッケージであったことが決め手 さらに解説 PDF を読むと Pickrell データも例題として使われていることが判明 3 をクリック 3 2 July 3,

69 GSVA をクリックすると 2Bioconductor から提供されているパッケージの場合はこんな感じになります 2 July 3,

70 Contents 機能解析 ( 発現変動遺伝子セット解析 ) 全体像 基本的な考え方と解析戦略の変遷 様々なプログラム 遺伝子セット情報の取得 (gmtファイルの取得) 発現データ情報と遺伝子セット情報のIDの対応付け 検証用 RNA-seqカウントデータセットPickrell data( なぜGSVAにしたか ) GSVAの解説 PDFを読み解く ( 手元のc.all.v6..entrez.gmt をどう読み込ませるか ) GSVAdataパッケージ提供の MSigDB c2コレクションであるc2broadsetsを理解する 手元のgmtファイルを読み込ませて GeneSetCollection 形式で取り扱えるようにする GSVA の解説 PDF を読み解く ( 手元の発現データファイルをどう取り扱うか ) ExpressionSetの取り扱い nsfilter 関数を用いた同一 IDの重複除去 メインプログラムgsva 関数が入力として受け付けるデータ形式 (ExpressionSetとMatrix) 検証用 RNA-seqカウントデータセットPickrell dataのイントロ スルーしていいところ MSigDB c2コレクションに2つの性特異的遺伝子セットを追加したものでgsvaを実行 ユニークなEntrez gene IDで グループごとに分離させた発現データファイル作成 整形後の発現データファイルとc.all.v6..entrez.gmtを入力としてGSVAを実行 July 3,

71 GSVA 左側の赤枠部分を拡大表示 July 3, 208 7

72 GSVA GSVA のインストール手順の説明 2Reference Manual の PDF これは CRAN( しーらん ) にもあるやつです 2 July 3,

73 GSVA の解説 PDF GSVA の解説 PDF(bignettes; ビニェット ) これは Bioconductor のサイトに存在するものです GSVA インストール後は R 起動直後に 2 をコピペでも ( 手元の PC 内に存在する ) 解説 PDF を開けます 2 July 3,

74 GSVA の解説 PDF GSVA の解説 PDF(bignettes; ビニェット ) これは Bioconductor のサイトに存在するものです GSVA インストール後は R 起動直後に 2 をコピペでも ( 手元の PC 内に存在する ) 解説 PDF を開けます 3 こんな感じ 3 2 July 3,

75 GSVA の解説 PDF GSVA の解説 PDF(bignettes; ビニェット ) を起動して 利用法を学んでいく 例えば 2 のところを押すと 3 のところに飛びます 基本的には解説 PDF 内を順番にざっと読んで このパッケージの使用法のノリに慣れる 2 3 July 3,

76 Contents 機能解析 ( 発現変動遺伝子セット解析 ) 全体像 基本的な考え方と解析戦略の変遷 様々なプログラム 遺伝子セット情報の取得 (gmtファイルの取得) 発現データ情報と遺伝子セット情報のIDの対応付け 検証用 RNA-seqカウントデータセットPickrell data( なぜGSVAにしたか ) GSVAの解説 PDFを読み解く ( 手元のc.all.v6..entrez.gmt をどう読み込ませるか ) GSVAdataパッケージ提供の MSigDB c2コレクションであるc2broadsetsを理解する 手元のgmtファイルを読み込ませて GeneSetCollection 形式で取り扱えるようにする GSVA の解説 PDF を読み解く ( 手元の発現データファイルをどう取り扱うか ) ExpressionSetの取り扱い nsfilter 関数を用いた同一 IDの重複除去 メインプログラムgsva 関数が入力として受け付けるデータ形式 (ExpressionSetとMatrix) 検証用 RNA-seqカウントデータセットPickrell dataのイントロ スルーしていいところ MSigDB c2コレクションに2つの性特異的遺伝子セットを追加したものでgsvaを実行 ユニークなEntrez gene IDで グループごとに分離させた発現データファイル作成 整形後の発現データファイルとc.all.v6..entrez.gmtを入力としてGSVAを実行 July 3,

77 MSigDB の C2 Application のあたりまでページ下部に移動 このパッケージは 遺伝子セット情報として 2MSigDB の C2 コレクションを利用していることがわかる 2 July 3,

78 MSigDB の C2 それは GSVAdata というパッケージに含まれる 2 c2broadsets という名前の 3GeneSetCollection オブジェクトだということがわかる 3 2 July 3,

79 MSigDB の C2 これをコピーして R コンソール画面上で コマンドのみペースト すれば確認できる July 3,

80 c2broadsets の確認 こんな感じでコピーして 2 コマンドのみペースト 2 July 3,

81 c2broadsets の確認 こんな感じになり GeneSetCollection という形式の c2broadsets を無事読み込めました July 3, 208 8

82 Tips: コマンドのみペースト コマンドのみペースト とすることで > などの余分な文字が誤ってコマンドとして認識されないようにすることができます July 3,

83 Contents 機能解析 ( 発現変動遺伝子セット解析 ) 全体像 基本的な考え方と解析戦略の変遷 様々なプログラム 遺伝子セット情報の取得 (gmtファイルの取得) 発現データ情報と遺伝子セット情報のIDの対応付け 検証用 RNA-seqカウントデータセットPickrell data( なぜGSVAにしたか ) GSVAの解説 PDFを読み解く ( 手元のc.all.v6..entrez.gmt をどう読み込ませるか ) GSVAdataパッケージ提供の MSigDB c2コレクションであるc2broadsetsを理解する 手元のgmtファイルを読み込ませて GeneSetCollection 形式で取り扱えるようにする GSVA の解説 PDF を読み解く ( 手元の発現データファイルをどう取り扱うか ) ExpressionSetの取り扱い nsfilter 関数を用いた同一 IDの重複除去 メインプログラムgsva 関数が入力として受け付けるデータ形式 (ExpressionSetとMatrix) 検証用 RNA-seqカウントデータセットPickrell dataのイントロ スルーしていいところ MSigDB c2コレクションに2つの性特異的遺伝子セットを追加したものでgsvaを実行 ユニークなEntrez gene IDで グループごとに分離させた発現データファイル作成 整形後の発現データファイルとc.all.v6..entrez.gmtを入力としてGSVAを実行 July 3,

84 c.all.v6..entrez.gmt 今手元にあるのは c.all.v6..entrez.gmt gmt ファイルを読み込ませ 遺伝子セット情報を自在に変えて解析できるようになりたい! July 3,

85 ?c2broadsets?c2broadsets と打ち込んで このデータがどのような手順で作成されたのか手がかりを探る July 3,

86 ?c2broadsets こんな感じになります GSVAdata パッケージ中の 2c2BroadSets の説明のページという意味 2 July 3,

87 ?c2broadsets c2broadsets という名前の 2GeneSetCollection 形式のオブジェクトは 3GSEABase パッケージ内の 4 getgmt 関数を用いて 5c2.all.v3.0.entrez.gmt ファイルを読み込んで得られたものだということが分かる July 3,

88 gmt ファイルの読込 2 例題 の gmt ファイル読込の基本形を実行してみましょう 2 July 3,

89 例題 : 基本形 コピペ実行結果 無事 GeneSetCollection 形式の geneset オブジェクトが得られていることがわかります 2 確かに C コレクションは 326 遺伝子セットでした 3 つめの遺伝子セットが chr5q23 で 42 つめが chr6q24 となっています July 3,

90 c.all.v6..symbols.gmt 入力ファイル (c.all.v6..symbols.gmt) を Excel で眺めたところ 2 最初の 2 行の遺伝子セット名と同じであり 妥当ですね 2 July 3,

91 c.all.v6..symbols.gmt これは 入力が c.all.v6..symbols.gmt なので 2PTMAP2 や 3FTMT のような gene symbols で遺伝子セット情報が記載されている 2 3 July 3, 208 9

92 c.all.v6..symbols.gmt c.all.v6..symbols.gmt 内にある gene symbols は 2PTMAP2や3FTMTを含めて全部で何種類あるのだろうか? そのあたりの情報は 2 3 July 3,

93 例題 : 基本形 c.all.v6..symbols.gmt 内には 2PTMAP2 や 3 FTMT を含めて 4 全部で 30,00 種類あるのだろう July 3,

94 まだ不十分か?! 赤下線部分に着目! この部分が NullIdentifier や NullCollection となっている GSVAdata パッケージの 2c2BroadSets ではそうなっていなかった 2 July 3,

95 両者を比較 gmt ファイルから読み込んだ geneset では 2NullIdentifier や NullCollection となっている その一方で GSVAdata パッケージの 3c2BroadSets では 4EntrezIdentifier や BroadCollection となっている ここまでやっておく必要性については今のところ不明ではあるが 念のためやったのが例題 July 3,

96 例題 3 2 例題 3 です 2 July 3,

97 例題 3 例題 3 の 2 入力は c.all.v6..entrez.gmt 2 July 3,

98 例題 3 getgmt 関数実行時に 2geneIdType と 3collectionType オプションを与えて Entrez gene ID であることや Broad institute 提供の C コレクションであることを明示しておけば 2 3 July 3,

99 例題 3 c.all.v6..entrez.gmt を読み込んで得られた 2geneset オブジェクトの中身が 3EntrezIdentifier や BroadCollection になります これで見た目上は GSVAdata パッケージの c2broadsets 同じような見栄えになりました 実際問題としてここまでやっておく必要があるかどうかはわかりません ここまでで gmt ファイルを読み込んで GeneSetCollection オブジェクトを作成するところまで完了 2 3 July 3,

100 Contents 機能解析 ( 発現変動遺伝子セット解析 ) 全体像 基本的な考え方と解析戦略の変遷 様々なプログラム 遺伝子セット情報の取得 (gmtファイルの取得) 発現データ情報と遺伝子セット情報のIDの対応付け 検証用 RNA-seqカウントデータセットPickrell data( なぜGSVAにしたか ) GSVAの解説 PDFを読み解く ( 手元のc.all.v6..entrez.gmt をどう読み込ませるか ) GSVAdataパッケージ提供の MSigDB c2コレクションであるc2broadsetsを理解する 手元のgmtファイルを読み込ませて GeneSetCollection 形式で取り扱えるようにする GSVA の解説 PDF を読み解く ( 手元の発現データファイルをどう取り扱うか ) ExpressionSetの取り扱い nsfilter 関数を用いた同一 IDの重複除去 メインプログラムgsva 関数が入力として受け付けるデータ形式 (ExpressionSetとMatrix) 検証用 RNA-seqカウントデータセットPickrell dataのイントロ スルーしていいところ MSigDB c2コレクションに2つの性特異的遺伝子セットを追加したものでgsvaを実行 ユニークなEntrez gene IDで グループごとに分離させた発現データファイル作成 整形後の発現データファイルとc.all.v6..entrez.gmtを入力としてGSVAを実行 July 3,

101 GSVA の解説 PDF 4. Functional enrichment のところ 最初に発現データとして 2 マイクロアレイデータの leukemia_eset を見せている これは 3ExpressionSet という発現データを格納する形式です 42,626 features 37 samples のデータのようですね July 3, 208 0

102 ExpressionSet GSVA パッケージでは RNA-seq カウントデータも ExpressionSet 形式になっています まだクリックしない! July 3,

103 重複除去時の入力 このあと行う同一 gene ID の重複除去時の入力として ExpressionSet 形式の 2leukemia_eset が与えられています 3 今は 6 ページのあたり 2 3 July 3,

104 重複除去時の入力 7~8 ページにかけて 2ExpressionSet オブジェクトの leukemia_eset を入力として 3nsFilter 関数を用いた重複除去が行われています 3 2 July 3,

105 nsfilter で重複除去 マイクロアレイ時代を知るヒトは AFFX という文字のみで 2leukemia_eset が Affymetrix GeneChip データであることがわかる また 3 の記述から Entrez gene ID であることを前提とし 4 で重複した Entrez gene ID の除去を行っているらしいことがわかる この段階で 重複除去を nsfilter 関数を用いて行うためには RNA-seq カウントデータの場合も ExpressionSet オブジェクトにしないといけないので 若干テンションが下がる July 3,

106 GSVA の入力 重複除去の実行結果は filtered_eset 28 ページ目の上のほうです 2 July 3,

107 Contents 機能解析 ( 発現変動遺伝子セット解析 ) 全体像 基本的な考え方と解析戦略の変遷 様々なプログラム 遺伝子セット情報の取得 (gmtファイルの取得) 発現データ情報と遺伝子セット情報のIDの対応付け 検証用 RNA-seqカウントデータセットPickrell data( なぜGSVAにしたか ) GSVAの解説 PDFを読み解く ( 手元のc.all.v6..entrez.gmt をどう読み込ませるか ) GSVAdataパッケージ提供の MSigDB c2コレクションであるc2broadsetsを理解する 手元のgmtファイルを読み込ませて GeneSetCollection 形式で取り扱えるようにする GSVA の解説 PDF を読み解く ( 手元の発現データファイルをどう取り扱うか ) ExpressionSetの取り扱い nsfilter 関数を用いた同一 IDの重複除去 メインプログラムgsva 関数が入力として受け付けるデータ形式 (ExpressionSetとMatrix) 検証用 RNA-seqカウントデータセットPickrell dataのイントロ スルーしていいところ MSigDB c2コレクションに2つの性特異的遺伝子セットを追加したものでgsvaを実行 ユニークなEntrez gene IDで グループごとに分離させた発現データファイル作成 整形後の発現データファイルとc.all.v6..entrez.gmtを入力としてGSVAを実行 July 3,

108 GSVA の入力 重複除去の実行結果は filtered_eset 28 ページ目の下のほうです 2 July 3,

109 ExpressionSet 形式 重複除去の実行結果は filtered_eset filtered_eset オブジェクト中の 2eset という部分の情報を抜き出した 3leukemia_filtered_eset が 4gsva 関数実行時の 5 入力のようです 3 leukemia_filtered_eset は ExpressionSet 形式です July 3,

110 GeneSetCollection 形式 c2broadsetsは GeneSetCollectionという形式の遺伝子セット情報です 2と3で解析する遺伝子セットのフィルタリングを指定しています 2は遺伝子セットを構成するメンバー数の下限 (minimum size) 3は上限 (maximum size) です どの遺伝子セット解析プログラムも 大抵このような遺伝子セットのフィルタリングを行います 従って 解析結果で見られる遺伝子セット数は 入力時よりも減るのが普通 2 3 July 3, 208 0

111 GSVA の入力形式 gsva の入力が 2ExpressionSet および 3GeneSetCollection という形式に限定されているかを 4?gsva で確認 July 3, 208

112 ?gsva こんな感じになります GSVA パッケージ中の 2gsva 関数の説明のページという意味 2 4 July 3, 208 2

113 ?gsva 何を書いてるのか (S4 method って何よ? とか ) 分かりづらいだろうが と 2 の比較から 2 4 July 3, 208 3

114 ?gsva 4 何を書いてるのか (S4 method って何よ? とか ) 分かりづらいだろうが と 2 の比較から 遺伝子セット情報は 3GeneSetCollection 形式以外に 4list 形式でもよいのだろう ということがわかる 3 2 July 3, 208 4

115 ?gsva 発現情報もまた ExpressionSet 以外に 2 matrix 形式でもよいことがわかる この結果から RNA-seq カウントデータの入力が通常のタブ区切りテキストファイルの場合は 基本そのまま読み込むのでよい ( 正確には as.matrix しないといけない ) と判断する 2 July 3, 208 5

116 おまけ 参考 遺伝子セットのフィルタリングは デフォルトでは行わない設定になっていることがわかる は遺伝子セットを構成するメンバー数の下限 ( minimum size) が 上限 (maximum size) が Inf になっているからです Inf は無限大の意味です July 3, 208 6

117 ?gsva 参考 このあたりにもちゃんと書いてますね July 3, 208 7

118 ?gsva kcdf オプションは この後の GSVA for RNA-seq data の記述を見てから気づくのが実際のところかもしれない 結論のみ述べると 2RNA-seq のカウントデータを入力とする場合は デフォルトの kcdf= Gaussian ではなく kcdf= Poisson で実行せねばならない 2 July 3, 208 8

119 Contents 機能解析 ( 発現変動遺伝子セット解析 ) 全体像 基本的な考え方と解析戦略の変遷 様々なプログラム 遺伝子セット情報の取得 (gmtファイルの取得) 発現データ情報と遺伝子セット情報のIDの対応付け 検証用 RNA-seqカウントデータセットPickrell data( なぜGSVAにしたか ) GSVAの解説 PDFを読み解く ( 手元のc.all.v6..entrez.gmt をどう読み込ませるか ) GSVAdataパッケージ提供の MSigDB c2コレクションであるc2broadsetsを理解する 手元のgmtファイルを読み込ませて GeneSetCollection 形式で取り扱えるようにする GSVA の解説 PDF を読み解く ( 手元の発現データファイルをどう取り扱うか ) ExpressionSetの取り扱い nsfilter 関数を用いた同一 IDの重複除去 メインプログラムgsva 関数が入力として受け付けるデータ形式 (ExpressionSetとMatrix) 検証用 RNA-seqカウントデータセットPickrell dataのイントロ スルーしていいところ MSigDB c2コレクションに2つの性特異的遺伝子セットを追加したものでgsvaを実行 ユニークなEntrez gene IDで グループごとに分離させた発現データファイル作成 整形後の発現データファイルとc.all.v6..entrez.gmtを入力としてGSVAを実行 July 3, 208 9

120 おさらい 2GSVA の解説 PDF には Pickrell データも例題として使われています 実際に Pickrell データを眺めます 3 2 July 3,

121 GSVA for RNA-seq data をクリックすると 2GSVA for RNA-seq data のところに飛びます 2 July 3, 208 2

122 p4 この画面内で実際に行うのは GSVAdata パッケージをロードした後 2 をコピペ実行することのみ 3 は無視でよい 2 3 July 3,

123 補足説明 の論文ではマイクロアレイデータが そして 2 の論文では RNAseq データが取得されており 両者は比較可能な状態にあります そして 2 の RNA-seq データはさらに Argonne sequencing center と Yale sequencing center の 2 か所で独立に取得されています 3 は単純に 4 アレイデータと ( 画面上では見えていませんが )p6 の 行目で見られる Argonne sequencing center で得られた RNA-seq データの gene ID が完全に一致しているかどうかを featurenames 関数で gene ID 情報を取り出した後 identical 関数で比較しているだけです 若干説明が不十分かもしれませんし私も誤解している部分があるかもしれませんが このあたりは深入りする価値はありません July 3,

124 p5~6 参考 のあたりが p6 の最初のほうの記述内容です これもさらっと流すところ July 3,

125 Contents 機能解析 ( 発現変動遺伝子セット解析 ) 全体像 基本的な考え方と解析戦略の変遷 様々なプログラム 遺伝子セット情報の取得 (gmtファイルの取得) 発現データ情報と遺伝子セット情報のIDの対応付け 検証用 RNA-seqカウントデータセットPickrell data( なぜGSVAにしたか ) GSVAの解説 PDFを読み解く ( 手元のc.all.v6..entrez.gmt をどう読み込ませるか ) GSVAdataパッケージ提供の MSigDB c2コレクションであるc2broadsetsを理解する 手元のgmtファイルを読み込ませて GeneSetCollection 形式で取り扱えるようにする GSVA の解説 PDF を読み解く ( 手元の発現データファイルをどう取り扱うか ) ExpressionSetの取り扱い nsfilter 関数を用いた同一 IDの重複除去 メインプログラムgsva 関数が入力として受け付けるデータ形式 (ExpressionSetとMatrix) 検証用 RNA-seqカウントデータセットPickrell dataのイントロ スルーしていいところ MSigDB c2コレクションに2つの性特異的遺伝子セットを追加したものでgsvaを実行 ユニークなEntrez gene IDで グループごとに分離させた発現データファイル作成 整形後の発現データファイルとc.all.v6..entrez.gmtを入力としてGSVAを実行 July 3,

126 p6 の上のほう 参考 のあたりも c2broadsets の中から 一部を抜き取っているだけ 個人的には なぜこれをやる必要があるのか理解不能 July 3,

127 C2 ではなく C 参考 c2broadsets の中から一部を抜き取って作成した canonicalc2broadsets は あくまでも C2 コレクションのデータの一部 その一方で 2 発現データの Pickrell data と MSigDB C コレクションの遺伝子セット解析で有意な発現変動を示したのは chryq であった July 3, AbsFilterGSEA(Yoon et al., PLoS One, : e06599, 206) 27

128 参考 性特異的セットを追加 canonicalc2broadsets のセットだけでは Pickrell data の解析結果を示しづらいので 22 つの sex-specific expression を示す遺伝子セットを追加すると書いています 2 July 3,

129 参考 性特異的セット つめ つめの sex-specific expression を示す遺伝子セットを MSY というセット名で作成したところ July 3,

130 参考 性特異的セット 2 つめ つめの sex-specific expression を示す遺伝子セットを MSY というセット名で作成したところ 22 つめの sex-specific expression を示す遺伝子セットを XiE というセット名で作成したところ 2 July 3,

131 参考 性特異的セットマージ 22つめのsex-specific expressionを示す遺伝子セットをxieというセット名で作成したところ 3 既存のcanonicalC2BroadSetsに対して 作成した2 つのsex-specific expressionを示す4msyと5 XiEを 6GeneSetCollectionという形式でマージした結果を 7 新たなcanonicalC2BroadSetsとする July 3, 208 3

132 Contents 機能解析 ( 発現変動遺伝子セット解析 ) 全体像 基本的な考え方と解析戦略の変遷 様々なプログラム 遺伝子セット情報の取得 (gmtファイルの取得) 発現データ情報と遺伝子セット情報のIDの対応付け 検証用 RNA-seqカウントデータセットPickrell data( なぜGSVAにしたか ) GSVAの解説 PDFを読み解く ( 手元のc.all.v6..entrez.gmt をどう読み込ませるか ) GSVAdataパッケージ提供の MSigDB c2コレクションであるc2broadsetsを理解する 手元のgmtファイルを読み込ませて GeneSetCollection 形式で取り扱えるようにする GSVA の解説 PDF を読み解く ( 手元の発現データファイルをどう取り扱うか ) ExpressionSetの取り扱い nsfilter 関数を用いた同一 IDの重複除去 メインプログラムgsva 関数が入力として受け付けるデータ形式 (ExpressionSetとMatrix) 検証用 RNA-seqカウントデータセットPickrell dataのイントロ スルーしていいところ MSigDB c2コレクションに2つの性特異的遺伝子セットを追加したものでgsvaを実行 ユニークなEntrez gene IDで グループごとに分離させた発現データファイル作成 整形後の発現データファイルとc.all.v6..entrez.gmtを入力としてGSVAを実行 July 3,

133 GSVA の実行準備完了 やっと GSVA を実行できる状況になったようです July 3,

134 p7 の上のほう が GSVA 実行コマンドです July 3,

135 p7 の上のほう Argonne sequence center で取得された RNAseq カウントデータの ExpressionSet オブジェクトと 2sex-specific expression を示す 2 つの遺伝子セット (MSY と XiE) を追加した MSigDB C2 コレクションの GeneSetCollection オブジェクトが入力 2 July 3,

136 p7 の上のほう 3 GSVA 実行結果である esrnaseq は 806 行 36 列からなるデータ 2 この後は RNA-seq カウントデータを RPKM 値に変換してから 3 マイクロアレイデータの結果である esmicro との比較を行って ほら似た結果になってるでしょ でオシマイ sexspecific expression を示す 2 つの遺伝子セット (MSY と XiE) については 一応サンプルごとに算出した Enrichment Scores (GSVA scores) の散布図を示している そして男女間でサンプルごとのスコアが確かに異なっており その傾向はマイクロアレイデータでも RPKM データでも同じですね ということは述べられている しかしながら 有意な発現変動遺伝子セットはどれかを調べる枠組みはガイドラインも示されておらず残念 2 July 3,

137 p7 の中央あたり RPKM 値を算出する最初のほうを示しているところ のあたりで RPKM 値への変換に必要な配列長情報を含む annoentrez22022 を取得 2 がまず RPM(Reads Per Million) 値を算出しているところ cpm 関数を用いていますが これは (Counts Per Million) 値を算出するものであり 実質的に同じです 3 は 4cpm と annotentrez22022 で共通の gene ID のもの (intersection) を抽出しているだけです July 3,

138 Contents 機能解析 ( 発現変動遺伝子セット解析 ) 全体像 基本的な考え方と解析戦略の変遷 様々なプログラム 遺伝子セット情報の取得 (gmtファイルの取得) 発現データ情報と遺伝子セット情報のIDの対応付け 検証用 RNA-seqカウントデータセットPickrell data( なぜGSVAにしたか ) GSVAの解説 PDFを読み解く ( 手元のc.all.v6..entrez.gmt をどう読み込ませるか ) GSVAdataパッケージ提供の MSigDB c2コレクションであるc2broadsetsを理解する 手元のgmtファイルを読み込ませて GeneSetCollection 形式で取り扱えるようにする GSVA の解説 PDF を読み解く ( 手元の発現データファイルをどう取り扱うか ) ExpressionSetの取り扱い nsfilter 関数を用いた同一 IDの重複除去 メインプログラムgsva 関数が入力として受け付けるデータ形式 (ExpressionSetとMatrix) 検証用 RNA-seqカウントデータセットPickrell dataのイントロ スルーしていいところ MSigDB c2コレクションに2つの性特異的遺伝子セットを追加したものでgsvaを実行 ユニークなEntrez gene IDで グループごとに分離させた発現データファイル作成 整形後の発現データファイルとc.all.v6..entrez.gmtを入力としてGSVAを実行 July 3,

139 発現行列データを整形して取得 順を追って説明をしますが 2 例題 7を実行しましょう 2 July 3,

140 例題 7 実行の下準備 話がややこしくなるので 一旦 R を再起動し 作業ディレクトリを Desktop に変更してから ls() を実行 これは利用可能なオブジェクトが何もないことを確認しているだけです 単に全員の環境を揃えているだけなので 意味が分かるヒトはやらなくてもよい July 3,

141 例題 7 赤枠内をコピペ GSVAdata パッケージ中にある発現データ commonpickrellhuang をロードし 2 どのようなオブジェクトが利用可能かを見るところまでです commonpickrellhuang には 3 種類の発現データが含まれています 3 アレイデータ および 4Argonne sequencing center と 5Yale sequencing center で取得された RNA-seq カウントデータです July 3, 208 4

142 列の並びはバラバラ 赤枠内をコピペ July 3,

143 列の並びはバラバラ Yale sequencing center で取得された RNA-seq カウントデータのオブジェクトを eset という名前で取り扱い 2 発現行列の列 (phenotype という意味で pdata) に関するメタデータ情報を確認したところ Gender 列が各サンプルの性別情報を含んでいると学習し そこを表示させた結果が 3 赤枠の最後のコマンド 列の並びは性別できれいに分離されていないことがわかる この段階で GSVA 実行結果を見やすくするために 予めソートしておこうと考える 2 3 July 3,

144 女 23 人 男 3 人 table 関数で内訳を調査 Female が 23 人 Male が 3 人 G 群 23 サンプル vs. G2 群 3 サンプル の 2 群間比較データだと読み替えてもよい July 3,

145 コード下部に移動 一番下までカーソル移動して最後の部分まで表示 July 3,

146 重複 gene ID の除去 オリジナルの発現行列 exprs(eset) に対して その行名からなるベクトル rownames(exprs(eset)) の要素数と 行名をユニークにした後のベクトルの要素数を表示 このデータの場合は = 26 個分の重複があったことがわかる この重複を除去してユニークな gene ID にしておかないと 保存 ( hoge7.txt の作成 ) 自体はできるが 次にそれを ( 私のウェブページ上の通常のやり方で ) 読み込むときに 重複した row.names は許されない! と言われる July 3,

147 nsfilter 関数の実行 重複除去を行う nsfilter 関数を実行 2 var.filter=f オプションがついているが こうしておかないと重複除去のみにならないからです 3 思い描いた通りの結果になっています 実際問題としては重複した Entrez gene IDs があった場合にどのような取り扱い ( 平均しているのかどうかなど ) を行っているかまでは調べきれていません 2 3 July 3,

148 列名変更の必要性 の後処理として列名変更をしなかった場合は 出力ファイル中の列名が 2 の実行結果のようになってしまいます 今は男女間での発現変動遺伝子セット解析を行いたいので 列名は Female or Male で充分 2 July 3,

149 列名変更後 赤枠内コピペ実行後 確かに列名変更できていることがわかります ただ 以前 pdata(eset)$gender でも確認したように Female と Male が性別順に並んでいないことがわかります July 3,

150 列名でソート 赤枠内コピペ実行後 確かに列名でソートできていることがわかります order 関数を用いた行列のソートテクは よく利用します July 3,

151 例題 7 を最後までコピペ後 hoge7.txt が作成されます July 3, 208 5

152 ,482 gene IDs hoge7.txt これが Entrez gene ID の重複除去を行って 性別ごとに列をソートした後の発現行列データ G 群 :23 females G2 群 :3 males July 3,

153 Contents 機能解析 ( 発現変動遺伝子セット解析 ) 全体像 基本的な考え方と解析戦略の変遷 様々なプログラム 遺伝子セット情報の取得 (gmtファイルの取得) 発現データ情報と遺伝子セット情報のIDの対応付け 検証用 RNA-seqカウントデータセットPickrell data( なぜGSVAにしたか ) GSVAの解説 PDFを読み解く ( 手元のc.all.v6..entrez.gmt をどう読み込ませるか ) GSVAdataパッケージ提供の MSigDB c2コレクションであるc2broadsetsを理解する 手元のgmtファイルを読み込ませて GeneSetCollection 形式で取り扱えるようにする GSVA の解説 PDF を読み解く ( 手元の発現データファイルをどう取り扱うか ) ExpressionSetの取り扱い nsfilter 関数を用いた同一 IDの重複除去 メインプログラムgsva 関数が入力として受け付けるデータ形式 (ExpressionSetとMatrix) 検証用 RNA-seqカウントデータセットPickrell dataのイントロ スルーしていいところ MSigDB c2コレクションに2つの性特異的遺伝子セットを追加したものでgsvaを実行 ユニークなEntrez gene IDで グループごとに分離させた発現データファイル作成 整形後の発現データファイルとc.all.v6..entrez.gmtを入力としてGSVAを実行 July 3,

154 C コレクションで GSVA 2 遺伝子セット情報は MSigDB C コレクションの c.all.v6..entrez.gmt を利用します これは 326 遺伝子セットからなります 2 July 3,

155 C コレクションで GSVA 発現データファイルの中身は July 3,

156 ,482 gene IDs SRP00540_23_3.txt Entrez gene ID の重複除去を行い 性別ごとに列をソートした後の発現行列データ G 群 :23 females G2 群 :3 males July 3,

157 コピペ実行 コード全体をコピペ実行後 全部で 326 個の遺伝子セットの発現変動解析を行うべく 入力として与えたが 2 遺伝子セットのメンバー数が 5 以上 500 以下という条件でフィルタリングすると 3298 個の遺伝子セットになったようです 3 2 July 3,

158 実行結果ファイル 実行結果ファイルの hoge.txt をエクセルで眺めてみましょう July 3,

159 hoge.txt GSVA 自体は Enrichment score を返すだけのプログラム そのスコアからなる数値ベクトルを入力として ノンパラメトリックな Wilcoxon rank sum test (Mann-Whitney U test と同じもの ) で得られた p 値が最も低い発現変動遺伝子セットは 2chryq でした この結果は 2 July 3,

160 同じ結果が得られた AbsFilterGSEA の論文の結果と同じですね 2 chryq がおそらく最上位だと思われます 3 解析データがちょっと異なりますが それでも結果は同じ 今回の我々の実行結果も 確かに male で高発現になっています July 3, 208 AbsFilterGSEA(Yoon et al., PLoS One, : e06599, 206) 60

161 hoge.txt male で高発現というのは の部分の Enrichment score が male 群で高いという理解で正しいはずです July 3, 208 6

162 Contents 機能解析 ( 発現変動遺伝子セット解析 ) July 3, 208 全体像 基本的な考え方と解析戦略の変遷 様々なプログラム 遺伝子セット情報の取得 (gmtファイルの取得) 発現データ情報と遺伝子セット情報のIDの対応付け 検証用 RNA-seqカウントデータセットPickrell data( なぜGSVAにしたか ) GSVAの解説 PDFを読み解く ( 手元のc.all.v6..entrez.gmt をどう読み込ませるか ) GSVAdataパッケージ提供の MSigDB c2コレクションであるc2broadsetsを理解する 手元のgmtファイルを読み込ませて GeneSetCollection 形式で取り扱えるようにする GSVA の解説 PDF を読み解く ( 手元の発現データファイルをどう取り扱うか ) ExpressionSetの取り扱い nsfilter 関数を用いた同一 IDの重複除去 メインプログラムgsva 関数が入力として受け付けるデータ形式 (ExpressionSetとMatrix) 検証用 RNA-seqカウントデータセットPickrell dataのイントロ スルーしていいところ MSigDB c2コレクションに2つの性特異的遺伝子セットを追加したものでgsvaを実行 ユニークなEntrez gene IDで グループごとに分離させた発現データファイル作成 整形後の発現データファイルとc.all.v6..entrez.gmtを入力としてGSVAを実行 最後に 62

163 最後に rcode_pickrell.txt 内で Pickrell カウントデータファイル (SRP00540_23_3.txt) を入力として以下の三つを行っています : サンプル間クラスタリング ( 雄雌入り混じった状態 ) シルエットスコアの計算 ( 非常に 0 に近い値 ) DEG 検出 (0%FDR で 85 個を DEG と判定 ) このように一見 hopeless な印象を受ける結果であったとしても 遺伝子セット解析によってカナリ説得力のある結果が得られることもあるというのが面白いですね July 3,