DESeq2で非発現変動遺伝子を見つける

発現変動解析では多くの場合発現変動している遺伝子、いわゆるDEGを検出するのが目的です。 ですので、帰無仮説はA群の平均発現量=B群の平均発現量という帰無仮説を立てて統計検定を行います。

ところで、研究の目的によっては発現変動「していない」遺伝子、非発現変動遺伝子の方が重要といった場合もあるかもしれません。この場合、通常の統計検定で帰無仮説が棄却されなかった遺伝子を非発現変動遺伝子としてしまうのは誤りだと考えられます。そこで考えられるのが帰無仮説を変えるという方法です。

DESeq2では帰無仮説を柔軟に設定できるので、今回はこの機能を使ってみたいと思います。まずはコードを提示します。

library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData = cts,
                              colData = coldata,
                              design = ~ condition)
dds <- DESeq(dds)
res <- results(dds, lfcThreshold=.5, altHypothesis="lessAbs")
plotMA(res)

f:id:menseki:20200715204345p:plain

図で赤くなっているところは基本的に検定で帰無仮説が棄却された遺伝子を示しています。図の中央横線付近は2群間の発現の比が1(つまりその対数は0)であるところであり、大体そのあたりに色がついていることから、「非」発現変動遺伝子が得られていることがわかります。

コードで重要なのは6行目の results 関数の引数で、 lfcThreshold はlog fold change の閾値を設定しており、 altHypothesis は対立仮説の不等号の向きです。今回はそれぞれ 0.5"lessAbs" を指定しているので、対立仮説は次のようになります。

 | \beta | < 0.5

ここで  \beta は log fold change を表しています。つまり、log fold change が0.5より小さいと考えられる遺伝子を取ってくるための検定の設定になっていると言えます。log fold change はいわば発現量比の対数なので、0のとき、発現変動していないことを表し、この値の絶対値が大きいほど発現の差も大きいことになります。符号は変動の向きです。

このようにDESeq2で直接指定するのは対立仮説ですが、帰無仮説をうまく設定することで様々な発現量比の遺伝子を検定によって同定することができます。設定を変えることで片側検定のようなこともできるようです。詳しくはマニュアルに書いてありますのでそちらをご参照ください。