・・・・・¡

いろいろ

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

参考

edgeRで3群間比較してみた

edgeRでの多群間発現変動解析について、少し実験してみました。

今回はこのような仮想的なデータを用意しました。G1、G2、G3という3つのグループを含むという設定です。今回用いるデータはかなりノイズが少ない現実のデータとは少し遠いシミュレーションデータを用いています。ただ、発現変動解析の挙動を把握するにはちょうど良いと思います。

発現パターン 遺伝子数
1 G1 = G2 = G3 16000
2 G1 < G2 = G3 500
3 G1 > G2 = G3 500
4 G2 < G3 = G1 500
5 G2 > G3 = G1 500
6 G3 < G1 = G2 500
7 G3 > G1 = G2 500
8 G1 < G2 < G3 500
9 G3 < G2 < G1 500

1行目の"G1 = G2 = G3"は非発現変動遺伝子nonDEGを表します。また、2,3行目はG1対他で発現変動、 4,5行目はG2対他で発現変動、 6,7行目はG3対他で発現変動している遺伝子を表します。最後に、8,9行目は全てのグループ間で発現変動している遺伝子を表します。これらの発現変動のパターンを図に示すと以下の様になります。また、図で表現するといろいろ面倒なので、下図中の矢印で示すようにそれぞれテキストでの表現を対応させます。この表現の仕方は、baySeqというRパッケージに倣っています。

f:id:menseki:20191126145039p:plain

table(truth)
{G1,G2,G3} {G1},{G2,G3} {G2},{G3,G1} {G3},{G1,G2} {G1},{G2},{G3}
16000 1000 1000 1000 1000

下準備

今回はこの仮想的なデータを用いて、有名な発現変動解析用RパッケージedgeRを用いて3群間の比較を行っていきたいと思います。 まずは下準備として、DGEList関数によるedgeR用のオブジェクト生成と、calcNormFactors関数による正規化係数の計算を行います。

library(edgeR)
##> Loading required package: limma
group <- gl(3, 3, labels = c("G1", "G2", "G3"))
design <- model.matrix(~group)
y <- DGEList(counts = counts, group = gl(3, 3, labels = c("G1", "G2", "G3")))

y <- calcNormFactors(y)

続いてデータのばらつきの指標であるDispersionパラメータを推定します。

y <- estimateCommonDisp(y)
y <- estimateTagwiseDisp(y)

G1 対 G2 の exact test

まずは正確検定によるG1とG2の比較を行います。

et <- exactTest(y, pair = c("G1", "G2"))
res <- topTags(et, n = nrow(counts), sort.by = "none")
pred_g1g2 <- factor(ifelse(res@.Data[[1]]$FDR <= 0.05, "DEG", "nonDEG"))
table(pred_g1g2, truth)
{G1,G2,G3} {G1},{G2,G3} {G2},{G3,G1} {G3},{G1,G2} {G1},{G2},{G3}
DEG 115 741 744 4 786
nonDEG 15885 259 256 996 214

edgeRではp値とともにFDRという値が得られるのですが、一般的にこの値を用いて帰無仮説を棄却するかを判定していくことが多いと思うので、ここではFDRが0.05以下の遺伝子を発現変動遺伝子とします。その結果、G1とG2の間で変動がある遺伝子がDEGと判定されていることがわかります。表の見方はRのtable関数のヘルプを参照してください。DEGと判定されているものを○で示すと以下の図のようになります。

f:id:menseki:20191126145042p:plain

G2 対 G3 の exact test

同様にしてG2対G3の正確検定もやってみます。

et <- exactTest(y, pair = c("G2", "G3"))
res <- topTags(et, n = nrow(counts), sort.by = "none")
pred_g2g3 <- factor(ifelse(res@.Data[[1]]$FDR <= 0.05, "DEG", "nonDEG"))
table(pred_g2g3, truth)
{G1,G2,G3} {G1},{G2,G3} {G2},{G3,G1} {G3},{G1,G2} {G1},{G2},{G3}
DEG 125 1 725 781 784
nonDEG 15875 999 275 219 216

f:id:menseki:20191126145045p:plain

G3 対 G1 の exact test

さらに、G3対G1の正確検定もやってみます。

et <- exactTest(y, pair = c("G3", "G1"))
res <- topTags(et, n = nrow(counts), sort.by = "none")
pred_g3g1 <- factor(ifelse(res@.Data[[1]]$FDR <= 0.05, "DEG", "nonDEG"))
table(pred_g3g1, truth)
{G1,G2,G3} {G1},{G2,G3} {G2},{G3,G1} {G3},{G1,G2} {G1},{G2},{G3}
DEG 132 766 8 755 1000
nonDEG 15868 234 992 245 0

ここで{G1},{G2},{G3}の遺伝子が100%DEGと判定されています。これは、今回の仮想データに含まれる発現パターンのうち、{G1},{G2},{G3}に対応するものはG1<G2<G3またはG3<G2<G1であり、G1とG3が発現量の差が最も大きいからだと考えられます。

f:id:menseki:20191126145049p:plain

exact testの全ての結果を統合

以上の結果をまとめて表示すると以下のようになります。

plyr::adply(table(pred_g1g2, pred_g2g3, pred_g3g1, truth), 1:3)
pred_g1g2 pred_g2g3 pred_g3g1 {G1,G2,G3} {G1},{G2,G3} {G2},{G3,G1} {G3},{G1,G2} {G1},{G2},{G3}
DEG DEG DEG 0 0 2 1 608
nonDEG DEG DEG 12 0 4 677 176
DEG nonDEG DEG 16 654 2 3 178
nonDEG nonDEG DEG 104 112 0 74 38
DEG DEG nonDEG 14 1 628 0 0
nonDEG DEG nonDEG 99 0 91 103 0
DEG nonDEG nonDEG 85 86 112 0 0
nonDEG nonDEG nonDEG 15670 147 161 142 0

一番理解しやすいのは1行目の情報で、全ての比較でDEGと判定された遺伝子は{G1},{G2},{G3}を予測しているという結果を示しています。逆に最終行では全ての比較でnonDEGと表現された場合は{G1,G2,G3}を予測しているといえるでしょう。また、2行目のG1対G2ではnonDEGと予測され、その他の組み合わせではDEGと判定された遺伝子が{G3},{G1,G2}という発現パターンのDEGを予測しているという結果が示されています。3行目と5行目も同様です。残った4,6,7行目はDEGの判定の組み合わせが矛盾しているデータです。例えば4行目ではG1対G2とG2対G3両方で発現変動していないと判定されていてこの時点でG3とG1でも発現変動していないのが自然なのですが、G1とG3で発現変動していると判定しています。このことは、統計検定の非対称性に起因するのではないかというのが私の考えです。つまり、検定で帰無仮説が棄却されなかった、FDR>0.05のデータに関して、DEGではない(nonDEGである)といった積極的な結論は導けないということです。

GLMバージョンの結果も示しておきます。

GLMに基づくDispersionパラメータの推定とデータへの当てはめ

y <- estimateGLMCommonDisp(y,design)
y <- estimateGLMTrendedDisp(y,design)
y <- estimateGLMTagwiseDisp(y,design)
fit <- glmFit(y,design)

G1 対 G2 のLRT

lrt <- glmLRT(fit,coef=2)
res <- topTags(lrt, sort.by = "none", n = nrow(counts))
pred_coef2 <- factor(ifelse(res@.Data[[1]]$FDR <= 0.05, "DEG", "nonDEG"))
table(pred_coef2, truth)
{G1,G2,G3} {G1},{G2,G3} {G2},{G3,G1} {G3},{G1,G2} {G1},{G2},{G3}
DEG 131 753 761 6 791
nonDEG 15869 247 239 994 209

f:id:menseki:20191126145052p:plain

G2 対 G3 でLRT

lrt <- glmLRT(fit,contrast = c(0, -1, 1))
res <- topTags(lrt, sort.by = "none", n = nrow(counts))
pred_cont <- factor(ifelse(res@.Data[[1]]$FDR <= 0.05, "DEG", "nonDEG"))
table(pred_cont, truth)
{G1,G2,G3} {G1},{G2,G3} {G2},{G3,G1} {G3},{G1,G2} {G1},{G2},{G3}
DEG 144 1 742 791 795
nonDEG 15856 999 258 209 205

f:id:menseki:20191126145056p:plain

G3 対 G1 でLRT

lrt <- glmLRT(fit,coef=3)
res <- topTags(lrt, sort.by = "none", n = nrow(counts))
pred_coef3 <- factor(ifelse(res@.Data[[1]]$FDR <= 0.05, "DEG", "nonDEG"))
table(pred_coef3, truth)
{G1,G2,G3} {G1},{G2,G3} {G2},{G3,G1} {G3},{G1,G2} {G1},{G2},{G3}
DEG 145 771 12 764 1000
nonDEG 15855 229 988 236 0

f:id:menseki:20191126145059p:plain

ANOVA的な方法

lrt <- glmLRT(fit,coef=2:3)
res <- topTags(lrt, sort.by = "none", n = nrow(counts))
pred_coef23 <- factor(ifelse(res@.Data[[1]]$FDR <= 0.05, "DEG", "nonDEG"))
table(pred_coef23, truth)
{G1,G2,G3} {G1},{G2,G3} {G2},{G3,G1} {G3},{G1,G2} {G1},{G2},{G3}
DEG 229 835 822 831 1000
nonDEG 15771 165 178 169 0

f:id:menseki:20191126145103p:plain

LRTの全ての結果を統合

plyr::adply(table(pred_coef2, pred_cont, pred_coef3, truth), 1:3)
pred_coef2 pred_cont pred_coef3 {G1,G2,G3} {G1},{G2,G3} {G2},{G3,G1} {G3},{G1,G2} {G1},{G2},{G3}
DEG DEG DEG 0 0 2 1 621
nonDEG DEG DEG 15 0 4 688 174
DEG nonDEG DEG 16 664 5 3 170
nonDEG nonDEG DEG 114 107 1 72 35
DEG DEG nonDEG 17 1 649 2 0
nonDEG DEG nonDEG 112 0 87 100 0
DEG nonDEG nonDEG 98 88 105 0 0
nonDEG nonDEG nonDEG 15628 140 147 134 0
plyr::adply(table(pred_coef2, pred_cont, pred_coef3, pred_coef23, truth), 1:4)
pred_coef2 pred_cont pred_coef3 pred_coef23 {G1,G2,G3} {G1},{G2,G3} {G2},{G3,G1} {G3},{G1,G2} {G1},{G2},{G3}
DEG DEG DEG DEG 0 0 2 1 621
nonDEG DEG DEG DEG 15 0 4 688 174
DEG nonDEG DEG DEG 16 664 5 3 170
nonDEG nonDEG DEG DEG 57 94 1 58 35
DEG DEG nonDEG DEG 17 1 649 2 0
nonDEG DEG nonDEG DEG 63 0 75 77 0
DEG nonDEG nonDEG DEG 56 73 84 0 0
nonDEG nonDEG nonDEG DEG 5 3 2 2 0
DEG DEG DEG nonDEG 0 0 0 0 0
nonDEG DEG DEG nonDEG 0 0 0 0 0
DEG nonDEG DEG nonDEG 0 0 0 0 0
nonDEG nonDEG DEG nonDEG 57 13 0 14 0
DEG DEG nonDEG nonDEG 0 0 0 0 0
nonDEG DEG nonDEG nonDEG 49 0 12 23 0
DEG nonDEG nonDEG nonDEG 42 15 21 0 0
nonDEG nonDEG nonDEG nonDEG 15623 137 145 132 0

参考

Rでタピオカミルクティー

plot.new()
plot.window(c(-10, 10), c(-10, 10), asp = 1)
polygon(c(-3, -4, 4, 3, -3), c(-8, 1, 1, -8, -8), col = "#ead7a4", border = F)
polygon(c(-3, -3 - 10 / 9, 3 + 10 / 9, 3, -3), c(-8, 2, 2, -8, -8), lwd = 2)
segments(c(0.8, 1.5), c(-7.6, -7.61), c(1.8, 2.5), c(4.4, 4.39), lwd = 2)
points(runif(30, -2.5, 2.5), runif(30, -7.5, -5), pch = 21, col = 0, bg = 1, cex = 2.5)

f:id:menseki:20191122010049p:plain

3グループのRNA-seqカウントデータのシミュレーション

3グループデータのモデル化

前回に続きましてRNA-seqの記事になります。まずは、3グループの場合のそれぞれのサンプル、遺伝子のカウントデータの平均パラメータを数式で簡単にモデル化していきたいと思います。こちらのページ(https://bi.biopapyrus.jp/rnaseq/analysis/de-analysis/3g-edger.html)にあるように、一般化線形モデルにおけるのリンク関数と線形予測子の考え方を用いて、ある遺伝子のサンプルごとの平均パラメータをベクトル \lambda で表すと、3つのグループの比較で各グループに3つの反復があるとき、以下のように考えることができます。

 
\log(\boldsymbol{\lambda}) = 
\begin{pmatrix}
1 & 0 & 0 \\
1 & 0 & 0 \\
1 & 0 & 0 \\
1 & 1 & 0 \\
1 & 1 & 0 \\
1 & 1 & 0 \\
1 & 0 & 1 \\
1 & 0 & 1 \\
1 & 0 & 1 \\
\end{pmatrix}
\begin{pmatrix}
\alpha \\
\beta_{1} \\
\beta_{2} \\
\end{pmatrix}

実際に計算して見ると\lambdaについて、以下のようになります。


\boldsymbol{\lambda} = 
\begin{pmatrix}
\exp(\alpha) \\
\exp(\alpha) \\
\exp(\alpha) \\
\exp(\alpha)\exp(\beta_1) \\
\exp(\alpha)\exp(\beta_1) \\
\exp(\alpha)\exp(\beta_1) \\
\exp(\alpha)\exp(\beta_2) \\
\exp(\alpha)\exp(\beta_2) \\
\exp(\alpha)\exp(\beta_2) \\
\end{pmatrix}

上から3つをGroup1、真ん中3つをGroup2、残り3つをGroup3のサンプルの平均とします。同じグループの反復は同じ平均を示していていい感じです。また、Group2はGroup1に対して\exp(\beta_1)倍だけ平均が異なり、Group3はGroup1に対して\exp(\beta_2)倍だけ異なる事がわかります。

2つのグループではGroup1とGroup2の発現の比(Fold change)を考えるだけでしたが、3つのグループを考えるときは2つのFold changeが考えられるということになります。

\alphaに関してはその遺伝子の全体的な発現レベルを決める役割を果たします。この値が小さければ全体的に発現が小さいということを表します。

では実際にこの式をもとに、前回の記事で考えた各発現パターンにおける平均パラメータ\lambdaを計算してみましょう。ここでは発現パターンについてのみ考えるので、\alpha = \log(100)で固定しています。

非発現変動遺伝子(nonDEG)

nonDEGの場合は簡単で、\beta_1 = 0, \beta_2 = 0とすれば良いです。

alpha <- log(100)
beta_1 <- log(1)
beta_2 <- log(1)
X <- model.matrix(~gl(3, 3))
lambda <- exp(X %*% c(alpha, beta_1, beta_2))
plot(lambda, type = "h", lwd = 30, lend = 1, col = "#696969", xlim = c(0, 10), ylim = c(-5, max(lambda)), axes = F, ann = F)
text(c(2, 5, 8), c(-3, -3, -3), c("Group1", "Group2", "Group3"), adj = c(0.5, 1))
segments(c(0.6, 3.6, 6.6), c(-1, -1, -1), c(3.4, 6.4, 9.4), c(-1, -1, -1))

f:id:menseki:20191121183611p:plain
nonDEG

発現変動遺伝子(1対他)

続いてあるグループでその他2つのグループと比べて発現変動している場合について考えます。Fold changeについては2とします。

G1 > G2 = G3

Group1においてのみ高発現の場合。

beta_1 <- log(1/2)
beta_2 <- log(1/2)
X <- model.matrix(~gl(3, 3))
lambda <- exp(X %*% c(alpha, beta_1, beta_2))
plot(lambda, type = "h", lwd = 30, lend = 1, col = "#696969", xlim = c(0, 10), ylim = c(-5, max(lambda)), axes = F, ann = F)
text(c(2, 5, 8), c(-3, -3, -3), c("Group1", "Group2", "Group3"), adj = c(0.5, 1))
segments(c(0.6, 3.6, 6.6), c(-1, -1, -1), c(3.4, 6.4, 9.4), c(-1, -1, -1))

f:id:menseki:20191121183653p:plain
G1 > G2 = G3

G1 < G2 = G3

Group1においてのみ低発現の場合。

beta_1 <- log(2)
beta_2 <- log(2)
X <- model.matrix(~gl(3, 3))
lambda <- exp(X %*% c(alpha, beta_1, beta_2))
plot(lambda, type = "h", lwd = 30, lend = 1, col = "#696969", xlim = c(0, 10), ylim = c(-5, max(lambda)), axes = F, ann = F)
text(c(2, 5, 8), c(-3, -3, -3), c("Group1", "Group2", "Group3"), adj = c(0.5, 1))
segments(c(0.6, 3.6, 6.6), c(-1, -1, -1), c(3.4, 6.4, 9.4), c(-1, -1, -1))

f:id:menseki:20191121221420p:plain
G1 < G2 = G3

G2 > G1 = G3

Group2においてのみ高発現の場合。

beta_1 <- log(2)
beta_2 <- log(1)
X <- model.matrix(~gl(3, 3))
lambda <- exp(X %*% c(alpha, beta_1, beta_2))
plot(lambda, type = "h", lwd = 30, lend = 1, col = "#696969", xlim = c(0, 10), ylim = c(-5, max(lambda)), axes = F, ann = F)
text(c(2, 5, 8), c(-3, -3, -3), c("Group1", "Group2", "Group3"), adj = c(0.5, 1))
segments(c(0.6, 3.6, 6.6), c(-1, -1, -1), c(3.4, 6.4, 9.4), c(-1, -1, -1))

f:id:menseki:20191121221520p:plain
G2 > G1 = G3

G2 < G1 = G3

Group2においてのみ低発現の場合。

beta_1 <- log(1/2)
beta_2 <- log(1)
X <- model.matrix(~gl(3, 3))
lambda <- exp(X %*% c(alpha, beta_1, beta_2))
plot(lambda, type = "h", lwd = 30, lend = 1, col = "#696969", xlim = c(0, 10), ylim = c(-5, max(lambda)), axes = F, ann = F)
text(c(2, 5, 8), c(-3, -3, -3), c("Group1", "Group2", "Group3"), adj = c(0.5, 1))
segments(c(0.6, 3.6, 6.6), c(-1, -1, -1), c(3.4, 6.4, 9.4), c(-1, -1, -1))

f:id:menseki:20191121221551p:plain
G2 < G1 = G3

G3 > G1 = G2

Group3においてのみ高発現の場合。

beta_1 <- log(1)
beta_2 <- log(2)
X <- model.matrix(~gl(3, 3))
lambda <- exp(X %*% c(alpha, beta_1, beta_2))
plot(lambda, type = "h", lwd = 30, lend = 1, col = "#696969", xlim = c(0, 10), ylim = c(-5, max(lambda)), axes = F, ann = F)
text(c(2, 5, 8), c(-3, -3, -3), c("Group1", "Group2", "Group3"), adj = c(0.5, 1))
segments(c(0.6, 3.6, 6.6), c(-1, -1, -1), c(3.4, 6.4, 9.4), c(-1, -1, -1))

f:id:menseki:20191121221621p:plain
G3 > G1 = G2

G3 < G1 = G2

Group3においてのみ低発現の場合。

beta_1 <- log(1)
beta_2 <- log(1/2)
X <- model.matrix(~gl(3, 3))
lambda <- exp(X %*% c(alpha, beta_1, beta_2))
plot(lambda, type = "h", lwd = 30, lend = 1, col = "#696969", xlim = c(0, 10), ylim = c(-5, max(lambda)), axes = F, ann = F)
text(c(2, 5, 8), c(-3, -3, -3), c("Group1", "Group2", "Group3"), adj = c(0.5, 1))
segments(c(0.6, 3.6, 6.6), c(-1, -1, -1), c(3.4, 6.4, 9.4), c(-1, -1, -1))

f:id:menseki:20191121221730p:plain
G3 < G1 = G2

全てのグループ間で発現変動

さらに続けて、全てのグループの間で発現変動している遺伝子について考えます。この場合もFolc changeについては
2を考えますが、一番発現が小さいグループと一番発現が多いグループのFolc changeは2^2 = 4となります。

G1 < G2 < G3
beta_1 <- log(2)
beta_2 <- log(4)
X <- model.matrix(~gl(3, 3))
lambda <- exp(X %*% c(alpha, beta_1, beta_2))
plot(lambda, type = "h", lwd = 30, lend = 1, col = "#696969", xlim = c(0, 10), ylim = c(-5, max(lambda)), axes = F, ann = F)
text(c(2, 5, 8), c(-3, -3, -3), c("Group1", "Group2", "Group3"), adj = c(0.5, 1))
segments(c(0.6, 3.6, 6.6), c(-1, -1, -1), c(3.4, 6.4, 9.4), c(-1, -1, -1))

f:id:menseki:20191121221923p:plain
G1 < G2 < G3

G1 < G3 < G2
beta_1 <- log(4)
beta_2 <- log(2)
X <- model.matrix(~gl(3, 3))
lambda <- exp(X %*% c(alpha, beta_1, beta_2))
plot(lambda, type = "h", lwd = 30, lend = 1, col = "#696969", xlim = c(0, 10), ylim = c(-5, max(lambda)), axes = F, ann = F)
text(c(2, 5, 8), c(-3, -3, -3), c("Group1", "Group2", "Group3"), adj = c(0.5, 1))
segments(c(0.6, 3.6, 6.6), c(-1, -1, -1), c(3.4, 6.4, 9.4), c(-1, -1, -1))

f:id:menseki:20191121221950p:plain
G1 < G3 < G2

G2 < G1 < G3
beta_1 <- log(1/2)
beta_2 <- log(2)
X <- model.matrix(~gl(3, 3))
lambda <- exp(X %*% c(alpha, beta_1, beta_2))
plot(lambda, type = "h", lwd = 30, lend = 1, col = "#696969", xlim = c(0, 10), ylim = c(-5, max(lambda)), axes = F, ann = F)
text(c(2, 5, 8), c(-3, -3, -3), c("Group1", "Group2", "Group3"), adj = c(0.5, 1))
segments(c(0.6, 3.6, 6.6), c(-1, -1, -1), c(3.4, 6.4, 9.4), c(-1, -1, -1))

f:id:menseki:20191121222034p:plain
G2 < G1 < G3

G3 < G1 < G2
beta_1 <- log(2)
beta_2 <- log(1/2)
X <- model.matrix(~gl(3, 3))
lambda <- exp(X %*% c(alpha, beta_1, beta_2))
plot(lambda, type = "h", lwd = 30, lend = 1, col = "#696969", xlim = c(0, 10), ylim = c(-5, max(lambda)), axes = F, ann = F)
text(c(2, 5, 8), c(-3, -3, -3), c("Group1", "Group2", "Group3"), adj = c(0.5, 1))
segments(c(0.6, 3.6, 6.6), c(-1, -1, -1), c(3.4, 6.4, 9.4), c(-1, -1, -1))

f:id:menseki:20191121222103p:plain
G3 < G1 < G2

G2 < G3 < G1
beta_1 <- log(1/4)
beta_2 <- log(1/2)
X <- model.matrix(~gl(3, 3))
lambda <- exp(X %*% c(alpha, beta_1, beta_2))
plot(lambda, type = "h", lwd = 30, lend = 1, col = "#696969", xlim = c(0, 10), ylim = c(-5, max(lambda)), axes = F, ann = F)
text(c(2, 5, 8), c(-3, -3, -3), c("Group1", "Group2", "Group3"), adj = c(0.5, 1))
segments(c(0.6, 3.6, 6.6), c(-1, -1, -1), c(3.4, 6.4, 9.4), c(-1, -1, -1))

f:id:menseki:20191121222127p:plain
G2 < G3 < G1

G3 < G2 < G1
beta_1 <- log(1/2)
beta_2 <- log(1/4)
X <- model.matrix(~gl(3, 3))
lambda <- exp(X %*% c(alpha, beta_1, beta_2))
plot(lambda, type = "h", lwd = 30, lend = 1, col = "#696969", xlim = c(0, 10), ylim = c(-5, max(lambda)), axes = F, ann = F)
text(c(2, 5, 8), c(-3, -3, -3), c("Group1", "Group2", "Group3"), adj = c(0.5, 1))
segments(c(0.6, 3.6, 6.6), c(-1, -1, -1), c(3.4, 6.4, 9.4), c(-1, -1, -1))

f:id:menseki:20191121222152p:plain
G3 < G2 < G1

Fold changeをランダムに設定する

これまではFold change(\beta)を固定して考えてきましたが、現実のRNA-seqデータをシミュレーションすることが目標なので、このパラメータをランダムに設定していきたいと思います。気をつけるべきところは、DEGの場合にランダムに設定してたまたまnonDEGとほぼ同じような発現データの状況になったということを避けることでしょう。そのために、発現変動があるというグループ間には
最低でも1.5倍の差が生じるようにしました。

G1 > G2 = G3
alpha <- log(100)
beta_1 <- beta_2 <- runif(1, min = -log(4), max = -log(1.5))
X <- model.matrix(~gl(3, 3))
lambda <- exp(X %*% c(alpha, beta_1, beta_2))
plot(lambda, type = "h", lwd = 30, lend = 1, col = "#696969", xlim = c(0, 10), ylim = c(-5, max(lambda)), axes = F, ann = F)
text(c(2, 5, 8), c(-3, -3, -3), c("Group1", "Group2", "Group3"), adj = c(0.5, 1))
segments(c(0.6, 3.6, 6.6), c(-1, -1, -1), c(3.4, 6.4, 9.4), c(-1, -1, -1))

f:id:menseki:20191121222228p:plain
G1 > G2 = G3

G1 < G2 = G3
alpha <- log(100)
beta_1 <- beta_2 <- runif(1, min = log(1.5), max = log(4))
X <- model.matrix(~gl(3, 3))
lambda <- exp(X %*% c(alpha, beta_1, beta_2))
plot(lambda, type = "h", lwd = 30, lend = 1, col = "#696969", xlim = c(0, 10), ylim = c(-5, max(lambda)), axes = F, ann = F)
text(c(2, 5, 8), c(-3, -3, -3), c("Group1", "Group2", "Group3"), adj = c(0.5, 1))
segments(c(0.6, 3.6, 6.6), c(-1, -1, -1), c(3.4, 6.4, 9.4), c(-1, -1, -1))

f:id:menseki:20191121222314p:plain
G1 < G2 = G3

G2 > G1 = G3
alpha <- log(100)
beta_1 <- runif(1, min = log(1.5), max = log(4))
beta_2 <- log(1)
X <- model.matrix(~gl(3, 3))
lambda <- exp(X %*% c(alpha, beta_1, beta_2))
plot(lambda, type = "h", lwd = 30, lend = 1, col = "#696969", xlim = c(0, 10), ylim = c(-5, max(lambda)), axes = F, ann = F)
text(c(2, 5, 8), c(-3, -3, -3), c("Group1", "Group2", "Group3"), adj = c(0.5, 1))
segments(c(0.6, 3.6, 6.6), c(-1, -1, -1), c(3.4, 6.4, 9.4), c(-1, -1, -1))

f:id:menseki:20191121222345p:plain
G2 > G1 = G3

G2 < G1 = G3
alpha <- log(100)
beta_1 <- runif(1, min = -log(4), max = -log(1.5))
beta_2 <- log(1)
X <- model.matrix(~gl(3, 3))
lambda <- exp(X %*% c(alpha, beta_1, beta_2))
plot(lambda, type = "h", lwd = 30, lend = 1, col = "#696969", xlim = c(0, 10), ylim = c(-5, max(lambda)), axes = F, ann = F)
text(c(2, 5, 8), c(-3, -3, -3), c("Group1", "Group2", "Group3"), adj = c(0.5, 1))
segments(c(0.6, 3.6, 6.6), c(-1, -1, -1), c(3.4, 6.4, 9.4), c(-1, -1, -1))

f:id:menseki:20191121222409p:plain
G2 < G1 = G3

G3 > G1 = G2
alpha <- log(100)
beta_1 <- log(1)
beta_2 <- runif(1, min = log(1.5), max = log(4))
X <- model.matrix(~gl(3, 3))
lambda <- exp(X %*% c(alpha, beta_1, beta_2))
plot(lambda, type = "h", lwd = 30, lend = 1, col = "#696969", xlim = c(0, 10), ylim = c(-5, max(lambda)), axes = F, ann = F)
text(c(2, 5, 8), c(-3, -3, -3), c("Group1", "Group2", "Group3"), adj = c(0.5, 1))
segments(c(0.6, 3.6, 6.6), c(-1, -1, -1), c(3.4, 6.4, 9.4), c(-1, -1, -1))

f:id:menseki:20191121222436p:plain
G3 > G1 = G2

G3 < G1 = G2
alpha <- log(100)
beta_1 <- log(1)
beta_2 <- runif(1, min = -log(4), max = -log(1.5))
X <- model.matrix(~gl(3, 3))
lambda <- exp(X %*% c(alpha, beta_1, beta_2))
plot(lambda, type = "h", lwd = 30, lend = 1, col = "#696969", xlim = c(0, 10), ylim = c(-5, max(lambda)), axes = F, ann = F)
text(c(2, 5, 8), c(-3, -3, -3), c("Group1", "Group2", "Group3"), adj = c(0.5, 1))
segments(c(0.6, 3.6, 6.6), c(-1, -1, -1), c(3.4, 6.4, 9.4), c(-1, -1, -1))

f:id:menseki:20191121222458p:plain
G3 < G1 = G2

G1 < G2 < G3
beta_1 <- runif(1, min = log(1.5), max = log(4))
beta_2 <- runif(1, min = beta_1 + log(1.5), max = beta_1 + log(4))
X <- model.matrix(~gl(3, 3))
lambda <- exp(X %*% c(alpha, beta_1, beta_2))
plot(lambda, type = "h", lwd = 30, lend = 1, col = "#696969", xlim = c(0, 10), ylim = c(-5, max(lambda)), axes = F, ann = F)
text(c(2, 5, 8), c(-3, -3, -3), c("Group1", "Group2", "Group3"), adj = c(0.5, 1))
segments(c(0.6, 3.6, 6.6), c(-1, -1, -1), c(3.4, 6.4, 9.4), c(-1, -1, -1))

f:id:menseki:20191121222523p:plain
G1 < G2 < G3

G1 < G3 < G2
beta_2 <- runif(1, min = log(1.5), max = log(4))
beta_1 <- runif(1, min = beta_2 + log(1.5), max = beta_2 + log(4))
X <- model.matrix(~gl(3, 3))
lambda <- exp(X %*% c(alpha, beta_1, beta_2))
plot(lambda, type = "h", lwd = 30, lend = 1, col = "#696969", xlim = c(0, 10), ylim = c(-5, max(lambda)), axes = F, ann = F)
text(c(2, 5, 8), c(-3, -3, -3), c("Group1", "Group2", "Group3"), adj = c(0.5, 1))
segments(c(0.6, 3.6, 6.6), c(-1, -1, -1), c(3.4, 6.4, 9.4), c(-1, -1, -1))

f:id:menseki:20191121222548p:plain
G1 < G3 < G2

G2 < G1 < G3
beta_1 <- runif(1, min = -log(4), max = -log(1.5))
beta_2 <- runif(1, min = log(1.5), max = log(4))
X <- model.matrix(~gl(3, 3))
lambda <- exp(X %*% c(alpha, beta_1, beta_2))
plot(lambda, type = "h", lwd = 30, lend = 1, col = "#696969", xlim = c(0, 10), ylim = c(-5, max(lambda)), axes = F, ann = F)
text(c(2, 5, 8), c(-3, -3, -3), c("Group1", "Group2", "Group3"), adj = c(0.5, 1))
segments(c(0.6, 3.6, 6.6), c(-1, -1, -1), c(3.4, 6.4, 9.4), c(-1, -1, -1))

f:id:menseki:20191121222608p:plain
G2 < G1 < G3

G3 < G1 < G2
beta_1 <- runif(1, min = log(1.5), max = log(4))
beta_2 <- runif(1, min = -log(4), max = -log(1.5))
X <- model.matrix(~gl(3, 3))
lambda <- exp(X %*% c(alpha, beta_1, beta_2))
plot(lambda, type = "h", lwd = 30, lend = 1, col = "#696969", xlim = c(0, 10), ylim = c(-5, max(lambda)), axes = F, ann = F)
text(c(2, 5, 8), c(-3, -3, -3), c("Group1", "Group2", "Group3"), adj = c(0.5, 1))
segments(c(0.6, 3.6, 6.6), c(-1, -1, -1), c(3.4, 6.4, 9.4), c(-1, -1, -1))

f:id:menseki:20191121222640p:plain
G3 < G1 < G2

G2 < G3 < G1
beta_2 <- runif(1, min = -log(4), max = -log(1.5))
beta_1 <- beta_2 + runif(1, min = -log(4), max = -log(1.5))
X <- model.matrix(~gl(3, 3))
lambda <- exp(X %*% c(alpha, beta_1, beta_2))
plot(lambda, type = "h", lwd = 30, lend = 1, col = "#696969", xlim = c(0, 10), ylim = c(-5, max(lambda)), axes = F, ann = F)
text(c(2, 5, 8), c(-3, -3, -3), c("Group1", "Group2", "Group3"), adj = c(0.5, 1))
segments(c(0.6, 3.6, 6.6), c(-1, -1, -1), c(3.4, 6.4, 9.4), c(-1, -1, -1))

f:id:menseki:20191121222704p:plain
G2 < G3 < G1

G3 < G2 < G1
beta_1 <- runif(1, min = -log(4), max = -log(1.5))
beta_2 <- beta_1 + runif(1, min = -log(4), max = -log(1.5))
X <- model.matrix(~gl(3, 3))
lambda <- exp(X %*% c(alpha, beta_1, beta_2))
plot(lambda, type = "h", lwd = 30, lend = 1, col = "#696969", xlim = c(0, 10), ylim = c(-5, max(lambda)), axes = F, ann = F)
text(c(2, 5, 8), c(-3, -3, -3), c("Group1", "Group2", "Group3"), adj = c(0.5, 1))
segments(c(0.6, 3.6, 6.6), c(-1, -1, -1), c(3.4, 6.4, 9.4), c(-1, -1, -1))

f:id:menseki:20191121222735p:plain
G3 < G2 < G1

疑似RNA-seqカウントデータ(3グループ)の生成

最後に3グループRNA-seqカウントデータのシミュレーションデータを生成してみたいと思います。ここまで、各発現パターンについて細かく考えてきましたが、シミュレーションデータを生成するだけなら、適当に\alpha, \beta_1, \beta_2を定めてやって各グループのデータの平均パラメータ\lambdaを求めてからそのベクトルの要素を適切に並び替えれば良いと思います。

あとはその平均パラメータをもとに負の二項分布に基づくノイズを加えてあげれば擬似的カウントデータの完成です。

以下には汎用性にかけますが、RNA-seqカウントデータ(3グループ)の生成コードを示します。

generateCountMeanMatrix <- function(Ngene, z) {
  X <- model.matrix(~gl(3,1))
  beta_1 <- runif(Ngene, min = log(1.5), max = log(4))
  beta_2 <- beta_1 + runif(Ngene, min = log(1.5), log(4))
  alpha <- rnorm(Ngene, 4, 1)
  tmp <- t(X %*% rbind(alpha, beta_1, beta_2))
  
  logMeanMat <- matrix(NA, nrow = Ngene, ncol = 9)
  # nonDEG
  logMeanMat[z ==  0,       ] <- tmp[z == 0, 1]
  # DEG type 1
  logMeanMat[z ==  1, 1:3   ] <- tmp[z == 1, 1]
  logMeanMat[z ==  1, -(1:3)] <- tmp[z == 1, 2]
  
  # DEG type 2
  logMeanMat[z ==  2,    1:3] <- tmp[z == 2, 2]
  logMeanMat[z ==  2, -(1:3)] <- tmp[z == 2, 1]
  
  # DEG type 3
  logMeanMat[z ==  3,    4:6] <- tmp[z == 3, 1]
  logMeanMat[z ==  3, -(4:6)] <- tmp[z == 3, 2]
  
  # DEG type 4
  logMeanMat[z ==  4,    4:6] <- tmp[z == 4, 2]
  logMeanMat[z ==  4, -(4:6)] <- tmp[z == 4, 1]
  
  # DEG type 5
  logMeanMat[z ==  5,    7:9] <- tmp[z == 5, 1]
  logMeanMat[z ==  5, -(7:9)] <- tmp[z == 5, 2]
  
  # DEG type 6
  logMeanMat[z ==  6,    7:9] <- tmp[z == 6, 2]
  logMeanMat[z ==  6, -(7:9)] <- tmp[z == 6, 1]
  
  # DEG type 7
  logMeanMat[z ==  7,    1:3] <- tmp[z == 7, 1]
  logMeanMat[z ==  7,    4:6] <- tmp[z == 7, 2]
  logMeanMat[z ==  7,    7:9] <- tmp[z == 7, 3]
  
  # DEG type 8
  logMeanMat[z ==  8,    1:3] <- tmp[z == 8, 1]
  logMeanMat[z ==  8,    4:6] <- tmp[z == 8, 3]
  logMeanMat[z ==  8,    7:9] <- tmp[z == 8, 2]
  
  # DEG type 9
  logMeanMat[z ==  9,    1:3] <- tmp[z == 9, 2]
  logMeanMat[z ==  9,    4:6] <- tmp[z == 9, 1]
  logMeanMat[z ==  9,    7:9] <- tmp[z == 9, 3]
  
  # DEG type 10
  logMeanMat[z ==  10,   1:3] <- tmp[z == 10, 3]
  logMeanMat[z ==  10,   4:6] <- tmp[z == 10, 1]
  logMeanMat[z ==  10,   7:9] <- tmp[z == 10, 2]
  
  # DEG type 11
  logMeanMat[z ==  11,   1:3] <- tmp[z == 11, 2]
  logMeanMat[z ==  11,   4:6] <- tmp[z == 11, 3]
  logMeanMat[z ==  11,   7:9] <- tmp[z == 11, 1]
  
  # DEG type 12
  logMeanMat[z ==  12,   1:3] <- tmp[z == 12, 3]
  logMeanMat[z ==  12,   4:6] <- tmp[z == 12, 2]
  logMeanMat[z ==  12,   7:9] <- tmp[z == 12, 1]
  
  logFC <- logMeanMat[, c(1, 4, 7)] - logMeanMat[, 1]
  
  list(logFC = logFC, logMeanMat = logMeanMat)
}


x <- generateCountMeanMatrix(Ngene = 20000, z = rep(0:12, c(14000, rep(500, 12))))
counts <- matrix(rnbinom(20000 * 9, size = 10, mu = exp(x$logMeanMat)), 20000, 9)
head(counts)

##>      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
##> [1,]   45   47   52   39   52   28   68   33   41
##> [2,]   86   67   65   35   45   52   84   43   41
##> [3,]  106  197  140  200  101  196   66  126  161
##> [4,]    9    6   12   18   10   10    6   19   17
##> [5,]   24   38   30   63   46   53   27   46   25
##> [6,]  120  237  246  157  270  396  329  163  189

heatmap(counts, Rowv = NA, Colv = NA)

f:id:menseki:20191121222806p:plain

それっぽいデータが得られました。

RNA-seqで3つ以上のグループの比較

突然ですがRNA-seqのリードカウントデータによる発現解析についての記事です。

今回はサンプルのグループが3つ以上のときの解析方法について、少し考えたことを述べていこうと思います。

3つ以上を一般的に考えると頭がこんがらがるので今回は3つの場合と限定しておきます。

このときの各Feature(遺伝子やIsoform、extonとか)の発現の種類について考えていきます。

2つのグループのときは発現変動しているかいないか(発現変動遺伝子(DEG) or nonDEG)で分けられ、さらにDEGについてはその大小関係から、グループ1(G1)で高発現かグループ2(G2)で高発現かの2種類に分けられます。

ということでDEGを大小関係で分けるところまで考えてもたかだか3パターンで済みます。

では3つのグループの比較ではどうなるでしょうか。言葉で説明するのが面倒なので図で示してみます。

f:id:menseki:20191120233800p:plain
3つのグループを比較したときの発現パターン

灰色の線より上にある1段めはDEGであるかどうか、DEGならどのグループ間で発現変動していないか、どこでは発現変動しているかを考えたときのパターンです。計5つが考えられます。

さらに、それぞれのパターンで発現変動の大小関係も考慮すると灰色の線以下のパターンが考えられることになります。

発現変動解析では発現変動遺伝子を得ることができればあとの大小関係はそれぞれのグループの平均を比べればわかります。したがって、上図の灰色線より上5つのパターンについてそれぞれのFeatureがどのパターンに属するかを考えるのが解析のメインで難しいところになってくるのかなと思います。

ちなみに、4グループの場合は発現変動の有無を考えるだけで15種類のパターンが考えられます。

次回はこの発現パターンをもとに3つのグループのRNA-seqリードカウントデータの疑似データを生成してみたいと思います。

簡易的に遺伝子型を図示したい

生物系の研究室ではある系統の交雑後代の遺伝子型を複数のマーカーについて調べるという実験をすることがよく?あります。こういった実験ではエクセルなどの表に結果を記録していきます。行にサンプル、列にマーカーが並んでいて、各セルにそのマーカーの示す遺伝子座のそれぞれサンプルの遺伝子型を書き込みす。例えばこんな感じの表ができあがったりします。

サンプル名 マーカー1 マーカー2 マーカー3
サンプル1 A A A
サンプル2 A B H

ある系統AとBの交雑後代を見ているとして、系統Aアレル型のホモをAと系統Bアレル型のホモをB、それらのヘテロをHとして
表記しているイメージです。

これを図にすることをよく研究でやっていたのですが、いちいちパワポで作っていたのでとても面倒でした。ということでこの図の描画作業を簡単に自動化して公開してみました。研究室ではこの図のことをグラフィカルジェノタイプといっていたので、そのままの題名で公開しています。

公開先はこちら→https://sabe.shinyapps.io/GraphicalGenoType/

最初の画面はこんな感じ。

f:id:menseki:20191112020252p:plain
最初の画面

使い方としてはまずgenotype fileとして先程の形式の表をタブ区切りのファイルとして指定してあげます。すると以下の図のように表にある情報を色付けして表示します。

f:id:menseki:20191112020301p:plain
遺伝子型情報を含む表形式ファイルを指定したところ

さらに、distance fileとして各マーカー間の距離を記載したファイル(1行に1つの距離、計【マーカー数‐1】行)を指定してあげると目的の図が表示されるようになります。あとは図のサイズが調整でき、図の中に記載するマーカー間の距離の単位をcMもしくはbpから選ぶことができます。

f:id:menseki:20191112020427p:plain
最終的な結果

図の説明を簡単にするとそれぞれの横棒が染色体の一部をイメージしています。縦線はマーカーの座乗位置を示していて、その座乗位置の遺伝子型がどうなっているかを色で表しています。黒がA型ホモ、白がB型ホモ、そして灰色がヘテロです。

現状はあまりカスタマイズできず、利用可能場面もかなり限られているのでもう少しバージョンアップできればなあと思っています。

もし興味がありましたらぜひ使ってみてください。