RNA-seqで多群間解析する(前置き)
RNA-seqのカウントデータから、発現変動解析で多群間比較を行う場合、全ての組み合わせに対してペアワイズに検定したり、グループをいったんプールして比較する(見たことありませんが)といった工夫が考えられます。ただ、このような方法は検定が多重化したり、帰無仮説が棄却されなかった場合の対処についていろいろ考えないといけない気がします。
例えば、3群間比較で、G1、G2、G3の3つのグループを比較したいとします。ここでは目標として、G1とその他では発現変動しているけどG2とG3の間では発現変動していないという遺伝子を見つけたいとします。
まず考えられる方法としては、G1とG2、G2とG3、G3とG1という組み合わせで発現変動解析を行うことだと思います。そして、直感的には、G1とG2およびG3とG1での比較では帰無仮説が棄却され、発現変動遺伝子と結論づけたが、一方でG2とG3の比較では帰無仮説が棄却されなかったという遺伝子が先程の目標の遺伝子である気がします。
ここで、2つの疑問が生じます。1つは複数回の検定を行ったが、有意水準は補正しなくて良いのか?また、2つ目は帰無仮説が棄却されないときは帰無仮説を採択して(発現変動していないと結論づけて)良いのか?です。
私は統計に関して言えば素人ですが、この方法にはちょっと違和感を感じます。
この違和感をうまく解決しているように思えるRNA-seq発現変動解析用のRパッケージにbaySeq
というものがあります。(私がそう思っているだけなので、実際そうなのかは論文の理論を読んでみてください)
baySeqで発現変動解析(本題)
さて、本題ですが、baySeqでは仮説の確率が最終的な出力になります。 ここで、仮説というのは発現変動の有無です。
edgeR
やDESeq2
といった統計検定ベースの方法では、発現変動の有無で2つの仮説を立て、その仮説のもとでのデータが得られる確率(尤度)をもとに検定を行い、帰無仮説を仮定するにはあまりにも極端なデータであると考えられたとき、帰無仮説を棄却するといった方法だと思います。
一方でbaySeq
では仮説の事後確率、つまり、手元のデータが得られたときの仮説の確率を計算します。2群間比較では発現変動の有無で2つの仮説を立てますが、3群間比較以上の場合は発現変動の有無と発現変動しているグループの組合わせで5つの仮説が建てられます。
考えられる5つの仮説(発現変動の有無とその組み合わせ)は以下の図に示すような形になります。
それぞれの仮説のbaySeq
流の表記法を矢印の先に示しました。
今回は以下に示すように{G1,G2,G3}
を除く各仮説ごとに1000の遺伝子を含む、計20000の遺伝子のカウントデータからなるシミュレーションデータを用いてbaySeq
での発現変動解析の流れを追っていこうと思います。
発現パターン | 遺伝子数 | baySeqでの表記法 | 遺伝子数 | |
---|---|---|---|---|
1 | G1 = G2 = G3 | 16000 | {G1,G2,G3} | 16000 |
2 | G1 < G2 = G3 | 500 | {G1},{G2,G3} | 1000 |
3 | G1 > G2 = G3 | 500 | ||
4 | G2 < G3 = G1 | 500 | {G1,G3},{G2} | 1000 |
5 | G2 > G3 = G1 | 500 | ||
6 | G3 < G1 = G2 | 500 | {G1,G2},{G3} | 1000 |
7 | G3 > G1 = G2 | 500 | ||
8 | G1 < G2 < G3 | 167 | {G1},{G2},{G3} | 1000 |
9 | G1 < G3 < G2 | 167 | ||
10 | G2 < G1 < G3 | 167 | ||
11 | G2 < G3 < G1 | 167 | ||
12 | G3 < G1 < G2 | 166 | ||
13 | G3 < G2 < G1 | 166 |
入力データはcounts
という引数に格納しており、各遺伝子の真の仮説(発現パターン)はtruth
に格納しています。
# countsの中身 head(counts) ##> G1.1 G1.2 G1.3 G2.1 G2.2 G2.3 G3.1 G3.2 G3.3 ##> gene_00001 3 1 2 17 0 1 8 10 1 ##> gene_00002 25 147 27 31 55 30 81 224 45 ##> gene_00003 23 0 182 0 0 3 17 2 0 ##> gene_00004 41 7 20 8 7 16 2 21 34 ##> gene_00005 42 39 45 60 47 112 54 66 83 ##> gene_00006 129 113 106 167 79 184 114 25 74 tail(counts) ##> G1.1 G1.2 G1.3 G2.1 G2.2 G2.3 G3.1 G3.2 G3.3 ##> gene_19995 21 13 59 19 7 8 0 2 0 ##> gene_19996 127 85 181 52 46 57 14 8 15 ##> gene_19997 743 431 588 5 13 240 106 38 7 ##> gene_19998 395 72 66 53 23 107 48 18 23 ##> gene_19999 802 603 250 206 193 280 88 57 90 ##> gene_20000 42 50 82 9 26 38 9 27 6 str(counts) ##> num [1:20000, 1:9] 3 25 23 41 42 129 0 15 51 481 ... ##> - attr(*, "dimnames")=List of 2 ##> ..$ : chr [1:20000] "gene_00001" "gene_00002" "gene_00003" "gene_00004" ... ##> ..$ : chr [1:9] "G1.1" "G1.2" "G1.3" "G2.1" ... # truthの中身 str(truth) ##> Factor w/ 5 levels "{G1,G2,G3}","{G1},{G2,G3}",..: 1 1 1 1 1 1 1 1 1 1 ... ##> - attr(*, "names")= chr [1:20000] "gene_00001" "gene_00002" "gene_00003" "gene_00004" ... table(truth) ##> truth ##> {G1,G2,G3} {G1},{G2,G3} {G1,G3},{G2} {G1,G2},{G3} {G1},{G2},{G3} ##> 16000 1000 1000 1000 1000
準備
まず、パッケージのロード、baySeq
のオブジェクト生成など下準備をします。
# パッケージをロード suppressPackageStartupMessages(library(baySeq)) # 各サンプルがどのグループに所属するかを示すベクトルの作成 replicates <- gl(3, 3, labels = c("G1", "G2", "G3")) # baySeqでの解析の基本となるcountDataオブジェクトの生成 cd <- new("countData", data = counts, replicates = replicates) # 各遺伝子ごとのアノテーション情報の付加 # ここでは遺伝子IDを設定 cd@annotation <- data.frame(gene_id = rownames(counts), stringsAsFactors = FALSE)
仮説の設定
準備ができたら、仮説の設定を行います。baySeq
の論文やドキュメントではmodelといっているところの設定です。ここではallModels
関数を使って考えうる全てのモデルを設定していますが、手動で設定することもできます。事前にこの仮説(発現パターン、model)はデータに含まれないとわかっていれば、その仮説(発現パターン、model)を除いて解析したほうが良いです。このことは、除く仮説の事前確率が0であるとするのと同じイメージだと思います。詳細は公式のドキュメント等を確認してください。
cd <- allModels(cd) # 設定された仮説(model)は以下のように表示させられる groups(cd) ##> $`{G1,G2,G3}` ##> [1] 1 1 1 1 1 1 1 1 1 ##> Levels: 1 ##> ##> $`{G1,G2},{G3}` ##> [1] 1 1 1 1 1 1 2 2 2 ##> Levels: 1 2 ##> ##> $`{G1,G3},{G2}` ##> [1] 1 1 1 2 2 2 1 1 1 ##> Levels: 1 2 ##> ##> $`{G1},{G2,G3}` ##> [1] 1 1 1 2 2 2 2 2 2 ##> Levels: 1 2 ##> ##> $`{G1},{G2},{G3}` ##> [1] 1 1 1 2 2 2 3 3 3 ##> Levels: 1 2 3
メインの計算部分
以下の2つの関数を実行することで、メインの計算を行います。各関数のcl
という引数で並列実行することができます。かなり時間がかかるので、複数のCPUコアを搭載しているマシンでは並列処理すると良いと思います。
# おそらく各仮説に対応する統計モデルのパラメータの事前分布を得る cd <- getPriors.NB(cd, cl = NULL) ##> Warning in getPriors.NB(cd, cl = cl): no library sizes available; inferring libsizes using default settings ##> Finding priors...done. # おそらく得られた事前分布をもとに各仮説の事後確率を計算 cd <- getLikelihoods(cd, cl = NULL) ##> Finding posterior likelihoods...Length of priorReps:0 ##> Length of priorSubset:20000 ##> Length of subset:20000 ##> Length of postRows:20000 ##> Analysing part 1 of 1 ##> Preparing data.....................................................done. ##> Estimating likelihoods......done! ##> . ##> done. parallel::stopCluster(cl)
結果
topCountsは各仮説ごとの事後確率が高い順に結果を示してくれます。group
引数でどの仮説についての結果を示すかを指定します。
likes
は事後確率の列です(Likelihoodは尤度なので紛らわしいですが)。同時にFDR(False discovery rate)とFWER(Family-wise error rate)も得られます。
baySeq
で画期的なのは非発現変動{G1,G2,G3}
の事後確率も得られるということと、一回の解析で5種類全ての発現変動の事後確率を得ることができるということだと勝手に思っています。
# 今回の場合は2から10列目には入力データが入っているので[-(2:10)]で省略しています。 # (普段はつける必要ないです) topCounts(cd, group = "{G1,G2,G3}")[-(2:10)] ##> gene_id likes FDR.{G1,G2,G3} FWER.{G1,G2,G3} ##> 6326 gene_06326 0.9999743 2.573369e-05 2.573369e-05 ##> 11483 gene_11483 0.9999295 4.811220e-05 9.622260e-05 ##> 7268 gene_07268 0.9999045 6.392063e-05 1.917509e-04 ##> 13906 gene_13906 0.9999014 7.257870e-05 2.902849e-04 ##> 1275 gene_01275 0.9998559 8.688787e-05 4.343676e-04 ##> 15386 gene_15386 0.9998508 9.727325e-05 5.835029e-04 ##> 14671 gene_14671 0.9998424 1.058887e-04 7.409923e-04 ##> 3463 gene_03463 0.9998384 1.128506e-04 9.024565e-04 ##> 1523 gene_01523 0.9998259 1.196614e-04 1.076447e-03 ##> 14057 gene_14057 0.9998014 1.275592e-04 1.274873e-03 topCounts(cd, group = "{G1},{G2,G3}")[-(2:10)] ##> gene_id likes {G1},{G2,G3} FDR.{G1},{G2,G3} FWER.{G1},{G2,G3} ##> 16346 gene_16346 0.9929350 2>1 0.007065008 0.007065008 ##> 16984 gene_16984 0.9849844 1>2 0.011040310 0.021974534 ##> 16334 gene_16334 0.9834026 2>1 0.012892659 0.038207174 ##> 16498 gene_16498 0.9825726 2>1 0.014026336 0.054968687 ##> 16787 gene_16787 0.9805773 1>2 0.015105617 0.073323785 ##> 16505 gene_16505 0.9801723 1>2 0.015892625 0.091697615 ##> 16205 gene_16205 0.9758792 2>1 0.017068073 0.113606558 ##> 16032 gene_16032 0.9738923 2>1 0.018198022 0.136748222 ##> 16325 gene_16325 0.9713096 2>1 0.019363841 0.161515256 ##> 16538 gene_16538 0.9712818 1>2 0.020299276 0.185595017 topCounts(cd, group = "{G1,G3},{G2}")[-(2:10)] ##> gene_id likes {G1,G3},{G2} FDR.{G1,G3},{G2} FWER.{G1,G3},{G2} ##> 17731 gene_17731 0.9958040 2>1 0.004195965 0.004195965 ##> 17003 gene_17003 0.9879476 1>2 0.008124186 0.016197801 ##> 17004 gene_17004 0.9867033 1>2 0.009848346 0.029279089 ##> 19662 gene_19662 0.9855506 1>2 0.010998613 0.043305436 ##> 19219 gene_19219 0.9839470 2>1 0.012009485 0.058663232 ##> 17609 gene_17609 0.9806427 2>1 0.013234118 0.076884952 ##> 19755 gene_19755 0.9797467 2>1 0.014236853 0.095581046 ##> 17851 gene_17851 0.9795680 2>1 0.015011245 0.114060123 ##> 17825 gene_17825 0.9786210 2>1 0.015718770 0.133000607 ##> 17671 gene_17671 0.9782681 2>1 0.016320080 0.151842120 topCounts(cd, group = "{G1,G2},{G3}")[-(2:10)] ##> gene_id likes {G1,G2},{G3} FDR.{G1,G2},{G3} FWER.{G1,G2},{G3} ##> 18740 gene_18740 0.9917037 2>1 0.008296295 0.008296295 ##> 18955 gene_18955 0.9905153 2>1 0.008890489 0.017702291 ##> 18283 gene_18283 0.9904999 1>2 0.009093693 0.027034216 ##> 18034 gene_18034 0.9831615 1>2 0.011029887 0.043417471 ##> 18871 gene_18871 0.9814391 2>1 0.012536081 0.061172463 ##> 18814 gene_18814 0.9813870 2>1 0.013548893 0.078646815 ##> 18908 gene_18908 0.9802947 2>1 0.014428382 0.096802370 ##> 18366 gene_18366 0.9801756 1>2 0.015102886 0.114707738 ##> 6441 gene_06441 0.9793903 1>2 0.015714752 0.132953328 ##> 18832 gene_18832 0.9784544 2>1 0.016297840 0.151634398 topCounts(cd, group = "{G1},{G2},{G3}")[-(2:10)] ##> gene_id likes {G1},{G2},{G3} FDR.{G1},{G2},{G3} FWER.{G1},{G2},{G3} ##> 19016 gene_19016 0.9074748 3>2>1 0.09252521 0.09252521 ##> 19368 gene_19368 0.8092578 3>1>2 0.14163370 0.26561893 ##> 19754 gene_19754 0.8002016 2>1>3 0.16102193 0.41234710 ##> 19889 gene_19889 0.7936998 1>2>3 0.17234151 0.53358003 ##> 19585 gene_19585 0.7801122 1>3>2 0.18185076 0.63614007 ##> 19264 gene_19264 0.7790403 2>3>1 0.18836891 0.71653843 ##> 19856 gene_19856 0.7557093 1>2>3 0.19635774 0.78578547 ##> 19476 gene_19476 0.7268015 3>1>2 0.20596283 0.84430855 ##> 19805 gene_19805 0.7265426 2>1>3 0.21346223 0.88688353 ##> 19324 gene_19324 0.7116312 2>3>1 0.22095289 0.91950280
それでは、各仮説について、FDR <= 0.05となるような仮説を採択して正解と照らし合わせてみましょう。
# 各仮説においてFDR<=0.05となるような遺伝子IDを取得 genes_model1 <- topCounts(cd, group = "{G1,G2,G3}", FDR = 0.05)$gene_id genes_model2 <- topCounts(cd, group = "{G1},{G2,G3}", FDR = 0.05)$gene_id genes_model3 <- topCounts(cd, group = "{G1,G3},{G2}", FDR = 0.05)$gene_id genes_model4 <- topCounts(cd, group = "{G1,G2},{G3}", FDR = 0.05)$gene_id genes_model5 <- topCounts(cd, group = "{G1},{G2},{G3}", FDR = 0.05)$gene_id ##> Warning in .selectTags(cD, group, ordering, decreasing = decreasing, number = number, : No features were found using the cutoffs for likelihood, FDR and FWER specified # 一応重複がないか確認 any(duplicated(c(genes_model1, genes_model2, genes_model3, genes_model4, genes_model5))) ##> [1] FALSE # 予測をベクトルで表現する pred <- character(20000L) names(pred) <- names(truth) pred[genes_model1] <- "{G1,G2,G3}" pred[genes_model2] <- "{G1},{G2,G3}" pred[genes_model3] <- "{G1,G3},{G2}" pred[genes_model4] <- "{G1,G2},{G3}" pred[genes_model5] <- "{G1},{G2},{G3}" # どの仮説ともされなかった遺伝子については保留としておく pred[pred == ""] <- "on hold" pred <- factor(pred, levels = c("on hold", "{G1,G2,G3}", "{G1},{G2,G3}", "{G1,G3},{G2}", "{G1,G2},{G3}", "{G1},{G2},{G3}")) # 正解と予測の混同行列 table(pred, truth) ##> truth ##> pred {G1,G2,G3} {G1},{G2,G3} {G1,G3},{G2} {G1,G2},{G3} {G1},{G2},{G3} ##> on hold 3226 585 587 583 789 ##> {G1,G2,G3} 12765 368 367 366 182 ##> {G1},{G2,G3} 3 45 0 0 10 ##> {G1,G3},{G2} 2 0 46 0 11 ##> {G1,G2},{G3} 4 2 0 51 8 ##> {G1},{G2},{G3} 0 0 0 0 0
残念ながら{G1},{G2},{G3}
については予測が0になってしまいましたが、それ以外のものについては一応目的のものを予測しようとしていることがわかりました。
上記の様に、一般的な発現変動解析と同じくFDRによって発現変動遺伝子の同定を行うことができますが、せっかく各遺伝子について各仮説ごとの事後確率が得られているので、MAP推定的な方法(最大事後確率の仮説を採択)も試してみたいと思います。baySeq
ではcd@posteriors
に全ての事後確率の対数が格納されているので、これを使いたいと思います。
# 事後確率の対数 head(cd@posteriors) ##> {G1,G2,G3} {G1,G2},{G3} {G1,G3},{G2} {G1},{G2,G3} {G1},{G2},{G3} ##> [1,] -0.13793033 -4.050070 -3.628333 -2.521969 -5.390949 ##> [2,] -0.11318505 -3.954112 -2.529247 -5.093921 -6.221638 ##> [3,] -0.10385443 -4.213384 -3.065705 -3.349006 -6.162442 ##> [4,] -0.07059272 -4.536489 -3.111822 -4.450149 -6.683720 ##> [5,] -0.10592185 -4.775803 -4.412836 -2.544166 -6.561011 ##> [6,] -0.09382480 -2.882759 -4.418776 -3.902618 -6.615311 pred <- max.col(cd@posteriors) pred <- colnames(cd@posteriors)[pred] pred <- factor(pred, levels = c("{G1,G2,G3}", "{G1},{G2,G3}", "{G1,G3},{G2}", "{G1,G2},{G3}", "{G1},{G2},{G3}")) # 正解と予測の混同行列 table(pred, truth) ##> truth ##> pred {G1,G2,G3} {G1},{G2,G3} {G1,G3},{G2} {G1,G2},{G3} {G1},{G2},{G3} ##> {G1,G2,G3} 15721 750 749 743 508 ##> {G1},{G2,G3} 93 232 13 13 153 ##> {G1,G3},{G2} 87 6 233 9 140 ##> {G1,G2},{G3} 97 12 4 234 133 ##> {G1},{G2},{G3} 2 0 1 1 66
何となく分類できているようですが、かなりFPも多いなという気がします。まあ、あくまで適当なシミュレーションデータの結果なので、本来の性能や、詳しい理論については直接論文を参照してください。