baySeqで多群間比較

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では仮説の確率が最終的な出力になります。 ここで、仮説というのは発現変動の有無です。

edgeRDESeq2といった統計検定ベースの方法では、発現変動の有無で2つの仮説を立て、その仮説のもとでのデータが得られる確率(尤度)をもとに検定を行い、帰無仮説を仮定するにはあまりにも極端なデータであると考えられたとき、帰無仮説を棄却するといった方法だと思います。

一方でbaySeqでは仮説の事後確率、つまり、手元のデータが得られたときの仮説の確率を計算します。2群間比較では発現変動の有無で2つの仮説を立てますが、3群間比較以上の場合は発現変動の有無と発現変動しているグループの組合わせで5つの仮説が建てられます。

考えられる5つの仮説(発現変動の有無とその組み合わせ)は以下の図に示すような形になります。 それぞれの仮説のbaySeq流の表記法を矢印の先に示しました。

G1 G2 G3 = = = G1 G2 G3 = G1 G2 G3 = G1 G2 G3 = G1 G2 G3 {G1,G2,G3} {G1},{G2,G3} {G1,G3},{G2} {G1,G2},{G3} {G1},{G2},{G3}

今回は以下に示すように{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も多いなという気がします。まあ、あくまで適当なシミュレーションデータの結果なので、本来の性能や、詳しい理論については直接論文を参照してください

参考