機能ゲノム学（第6回）

Size: px

Start display at page:

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

そよかみいしづ
6 years ago
Views:

1 RNA-Seqデータ解析における正規化法の選択 :RPKM 値でサンプル間比較は危険?! 東京大学大学院農学生命科学研究科アグリバイオインフォマティクス教育研究ユニット門田幸二 ( かどたこうじ ) kadota@iu.a.u-tokyo.ac.jp 1

2 よりよい正規化法とは? その正規化法によって得られたデータを用いて発現変動の度合いでランキングしたときに真の発現変動遺伝子 ( DEG) がより上位にランキングされる ( 感度特異度高い ) 正規化法 1 正規化法 % 83.3% Area Under the ROC Curve (ROC 曲線の下部面積 :AUC) 2

3 Marioni et al., Genome Res., 18: , 2008 のデータ出発点 ( マイクロアレイと同じく ) マップされたリード数情報を含む遺伝子発現行列データ 3

代表的な正規化法 :RPKM RPM 正規化 ( マイクロアレイなどと同じところ ) Reads per million mapped reads サンプルごとにマップされた総リード ( 塩基配列 ) 数が異なる各遺伝子のマップされたリード数を総 read 数が100 万 (one million) だった場合に補正 raw counts:all reads= RPM :

4 代表的な正規化法 :RPKM RPM 正規化 ( マイクロアレイなどと同じところ ) Reads per million mapped reads サンプルごとにマップされた総リード ( 塩基配列 ) 数が異なる各遺伝子のマップされたリード数を総 read 数が100 万 (one million) だった場合に補正 raw counts:all reads= RPM : 1,000,000 A1BGの場合は 744 : 5,087,097 = RPM : 1,000,000 1,000,000 1,000,000 RPM raw counts all reads 5,087,097 RPKM 正規化 (RNA-seq 特有 ) Mortazavi et al., Nat. Methods, 5: , 2008 Reads per kilobase of exon per million mapped reads 遺伝子の配列長が長いほど配列決定 (sequence) される確率が上昇各遺伝子の配列長を 1000 塩基 (one kilobase) の長さだった場合に補正 RPKM 1,000,000 1,000 raw counts all reads gene length 1,000,000,000 A1BG ,764 5,087,097 1,000,000,000 raw counts gene length all reads 4

5 RPKM( 正確には RPM) の問題点仮定全 4 遺伝子長さが同じ (gene length の項を無視できるので ) 遺伝子 4 だけが発現変動遺伝子 (DEG) サンプル A (all reads = 15) 遺伝子 1 遺伝子 2 遺伝子 3 遺伝子 4 Robinson and Oshlack, Genome Biol., 11:R25, 2010 定数 raw counts gene length all reads サンプル A (all reads = 30) 遺伝子 1 遺伝子 2 遺伝子 3 遺伝子 4 補正サンプルB (all reads = 30) 遺伝子 1 遺伝子 2 遺伝子 3 遺伝子 4 サンプルB (all reads = 30) 遺伝子 1 遺伝子 2 遺伝子 3 遺伝子 4 補正後の解析結果 :A で高発現が 3 個, B で高発現が 1 個 5

6 M TMM 正規化法 RPM 補正後のデータ Robinson and Oshlack, Genome Biol., 11:R25, 2010 TMM 補正後のサンプル B のデータ P DEG 2 (%) TMM = A P DEG 2 縦軸で上位下位合わせて P DEG % を Trim 残りのデータで M の weighted Mean(TMM) を計算 (%) 6

7 TMM 補正の有無で結論が異なることも得られた DEG セット中の割合 TMM 補正なし (Marioni et al., Genome Res., 2008) サンプルA(Kidney):78% サンプルB(Liver):22% Robinson and Oshlack, Genome Biol., 11:R25, 2010 TMM 補正あり (Robinson and Oshlack, 2010) サンプル A(Kidney):53% サンプル B(Liver):47% TMM 法で使用されているパラメータ ( 一部 ) log 2 (B/A) で発現変動順にランキングし全体で全遺伝子数の 60% 分を Trim (P DEG = 60%) その内訳はサンプル A 側とサンプル B 側で高発現なものを各 50% とする (P A = 50%) A 群 B 群 PDEG P A P DEG ( 100 PA ) Trim 後に残ったデータのみを用いて正規化係数を決定 7

8 参考 URL 8

9 TMM-normalized data の作成法 9

10 利用可能な R パッケージたち DEGseq (Wang et al., Bioinformatics, 26: , 2010) ポワソン分布 (variance = mean) を仮定しているためばらつきを過少評価 edger (Robinson et al., Bioinformatics, 26: , 2010) 正規化法 :TMM 法負の二項分布 (variance > mean) を仮定 mean のみのパラメータを用いて現実のばらつきを表現 DESeq (Anders and Huber, Genome Biol., 11: R106, 2010) 正規化法 :RLE 法 (relative log expression) edger のモデルをさらに拡張 ( しているらしい ) bayseq (Hardcastle and Kelly, BMC Bioinformatics, 11:422, 2010) 正規化法 :RPM ( たぶん ) 配列の長さ情報を与えるオプションがあるデータセット中に占める DEG の割合 (P DEG ) を一意に返す NBPSeq (Di et al., SAGMB, 10:24, 2011) アルゴリズムの詳細ついては不明です 10

11 TMM 法と門田法 TMM 法で用いられているパラメータ P DEG = 60%( 全遺伝子数の 60% 分を Trim した残りのデータで正規化係数を決定 ) P A = 50%( 発現変動遺伝子の内訳 : A 群 > B 群 = A 群 < B 群 ) 門田法 bayseq で自動推定された P DEG に相当しないデータを用いて正規化係数を決める 11

12 シミュレーションデータ Robinson and Oshlack, Genome Biol., 11:R25, 2010 TMM paper の Fig. 3 で用いられたものと同じ関数を使用ポアソン分布発現変動遺伝子 (DEG) の倍率変化 : 2 Number of common genes: 20,000 sample Aのみで発現している遺伝子数 : 2,200 0 sample Bのみで発現している遺伝子数 : 0 DEGの割合 (P DEG ): 5% DEG 中に占めるsampleA > Bの割合 (P A ): 80% 全遺伝子数が 20,000 個そのうち 5%(1,000 個 ) が DEG(P DEG = 5%) 80%(800 個 ) が samplea で高発現 (P A = 80%) 12

13 門田法のイメージ真実 P DEG = 5%, P A = 80% bayseq による推定値 P DEG = 3.4%, P A = 73.3% bayseq で DEG と判定されなかった ( 灰色部分に相当 ) データのみを用いて正規化係数を決定 13

14 評価基準 AUC 値 ( 高いほど感度特異度が高いことを意味する ) あるシミュレーション条件 (P DEG =5%, P A =80%) の結果想定される様々なシナリオに対する頑健性はどうか? 14

15 P DEG 5% 様々なシナリオ P A = 50% 60% 70% 80% 90% 100% 10% 20% 30% 15

16 P DEG 5% 門田法のイメージ P A = 50% 60% 70% 80% 90% 100% 10% 20% 30% 16

17 リアルデータ Hardcastle and Kelly, BMC Bioinformatics, 11:422, 2010 bayseq 論文の Fig. 5 のデータ公共データベース (GEO): GSE16959 Arabidopsis thaliana の leaf samples の塩基の small RNAs (srnas) two wild-type (WT) samples vs. two RDR6 (RNA-dependent RNA polymerase 6) knockout (KO) samples RDR6 は tasrnas (trans-acting srnas) 生成に必要であることが既知 70,619 unique small RNA sequences が Arabidopsis thaliana ゲノムにヒット tasrna loci のみに 100% マッチし発現が WT > KO となる 657 potentially true positives を同定 (P A = 100%) 70,619 行 4 列の srna 発現行列中に 657 個の真の発現変動 srnas を含むデータ 17

18 解析結果 (ROC 曲線の一部 ) Fig. 5 in bayseq paper Hardcastle and Kelly, BMC Bioinformatics, 11:422,

19 まとめ性能評価結果シミュレーション : bayseq > DESeq > edger > NBPSeq あるリアルデータ : bayseq edger > NBPSeq > DESeq 正規化法はよりよいものを使うべし RPM 法, TMM 法, 門田法などオリジナルの手順よりも門田法と組み合わせるとよりよい結果に計算環境 Windows 7, Intel Xeon 3.33GHz (2 processor), 96GB メモリ 19

20 謝辞共同研究者西山智明博士 ( 金沢大学 ) 清水謙多郎博士 ( 東京大学 ) グラント若手研究 (B)(H21-23 年度 ): マイクロアレイ解析の再現性感度特異度を飛躍的に向上させるデータ解析手法の開発 ( 代表 ) 新学術領域研究 ( 研究領域提案型 )(H22 年度 -): 非モデル生物におけるゲノム解析法の確立 ( 分担 ; 研究代表者 : 西山智明 ) 解析データ提供 Dr. Thomas J. Hardcastle(baySeq 著者 ) 20

21 その他 ( 苦労した点など ) bayseq 実行時のリサンプリング回数を当初 1000 回にしていた結果がころころ変わって困ったよくよくマニュアルを見ると推奨は 10,000 回だった bayseq 原著論文の精度はひょっとして DESeq と bayseq は TMM や門田法で正規化した後のデータをどうやって入力データとすればいいのか最終的にこの二つは同じやり方でできた edger もしばらく上記二つと同じやり方でやっていたが最近 ( 9 月 ) そのやり方ではだめだということに気付いた一つのシミュレーション条件を 100 回で 30 時間かかります 21

22 その他 ( 実は一番重要かも ) 読み込ませる入力データ 1: 通常 ( 生のリードカウント ) TMM 正規化係数などは引数で与える 2:TMM or 門田正規化後のデータプログラム内部で正規化されないように操作 3:RPM 正規化後のデータ TMM 正規化係数などは引数で与えるリアルデータ 22

23 シミュレーションデータ NBPSeq のデータを生データの (μ の ) 経験分布として使用負の二項分布 ( 平均 =μ 分散 =μ+μ 2 /size;size=10とした) 発現変動遺伝子 (DEG) の倍率変化 : 4 Number of common genes: 20,000 DEGの割合 (P DEG ): 5, 10, 20, 30% DEG 中に占めるsampleA > Bの割合 (P A ): 50, 60, 70, 80, 90, 100% 23

24 P DEG 5% 様々なシナリオ P A = 50% 60% 70% 80% 90% 100% 10% 20% 30% 24

25 P DEG 5% 様々なシナリオ P A = 50% 60% 70% 80% 90% 100% 10% 20% 30% 25

機能ゲノム学（第6回）

機能ゲノム学（第6回）トランスクリプトーム解析の今昔なぜマイクロアレイ? なぜRNA-Seq? 東京大学大学院農学生命科学研究科アグリバイオインフォマティクス教育研究ユニット門田幸二 ( かどたこうじ ) http://www.iu.a.u-tokyo.ac.jp/~kadota/ kadota@iu.a.u-tokyo.ac.jp 1 Contents トランスクリプトーム解析の概要各手法の長所短所マイクロアレイ