寒冷ストレス処理群と対照群の比較実験を行い、両者の間で発現量が異なる発現変動遺伝子について考える。実験デザインとして、ストレス処理群が 100 biological replicates、対照群が 100 biological replicates であるとする。ここで、正規化後のデータを利用して、gene 1 のストレス処理群と対照群それぞれにおける 100 replicates の平均発現量を計算し、グラフで可視化する。また、別の gene 2 についても同様に平均を計算して可視化する。可視化した結果を見ると、gene 1 と gene 2 の両方とも、ストレス処理群での発現量が 2 倍高いことがわかる。このとき、gene 1 と gene 2 のストレス処理群での発現量が、対照群での発現量に比べて有意に異なると言えるのだろうか。正確に言えば、この図だけでは判断できない。

ここで上の図にそれぞれの実験群の 100 replicates のデータを追加して可視化する。この 100 replicates の測定値からヒストグラムを描くことができ、さらにその発現量の分布（平均と分散）をほぼ正確に推定できるようになる。このとき、それぞれの遺伝子に対して、処理群と実験群における発現量の分布の重なりに基づいて、発現量の平均に有意差があるかどうかを調べることができる。

一つの実験群で 100 replicates を用意するのは非現実的である。そこで、一つの実験群で 3 replicates を用意して実験を行ったと考えてみる。replicate が少なくなると、各群の平均と分散を正確に求めることはできないが、ある程度、真に近い値を推定できる。このとき、gene 1 と gene 2 の発現量の分布（平均値と分散）を見比べてみる。gene 1 では、ストレス処理群と対照群の分布が互いに離れており、オーバーラップがほとんどない。そのため、再実験を繰り返したとしても、ストレス処理群であるにもかかわらず対照群と同じ発現量になるようなサンプルが現れることがない（95% 以上の確率で現れることがない）。この場合、gene 1 の発現量は、ストレス処理群と対照群の間で有意に差が見られると言える。これに対して、gene 2 の分布を確認すると、ストレス処理群の分布と対照群の分布が重なっている部分が大きい。このため、実験を繰り返えすと、ストレス処理群であるにもかかわらず対照群と同じ発現量になるようなサンプルが頻繁に現れるようになる。この場合、gene 2 がストレス処理群のサンプルと対照群のサンプルで同じ発現量をとるケースが多く観測できるので、両者の間で有意差はないと言える。

上で見たように発現変動遺伝子を検出するために、分散を正確に予測することが非常に重要である。この分散を予測するためには、replicates が必要である。したがって、コストが高くなるものの、RNA-Seq の実験では必ず replicates を入れべきである。replicates のないデータでは分散を予測できず、その検定結果も信用できない。現在では、発現変動遺伝子を検出するのが目的であれば、各実験群で 3 replicates を用意するのが一般的である。

分散の推定

RNA-Seq のリードカウントの平均 μ と分散 V を求めて、散布図で可視化すると下図のようになる。RNA-Seq データは、分散が平均よりも大きい。分散と平均が同じであれば、ポアソン分布でモデル化して、平均を使って分散を推定すればよいが、RNA-Seq データの場合は、分散と平均の間にそのような単純な関係がないことがわかる。このとき、分散が平均よりも大きくなるように補正を行えばよく、その補正を行う因子をここで定数 φ とおく。補正方法は様々な形で考えられるが、例えば、次のように補正して、平均と分散の関係を式で表してみることにする。

定数 φ を用いた補正では、平均 μ が小さいときに \(V = \mu + \mu^{2}\phi\) が正確に平均と分散の関係を表せるが、平均が大きくなると実際のデータとこの式の間にズレが生じるようになる。そこで、φ を定数ではなく、変数と考える。平均が小さいとき φ が大きな値をとり、平均が大きくなったときに φ が小さくなるような変数と考える。このとき、あらためて式 \(V = \mu + \mu^{2}\phi\) と実際のデータを可視化すると、以下のようにズレがなくなる。つまり、分散 V を正確に推定するためには、平均 μ とパラメーター φ の関係を明らかにしなければならない。この φ のことを dispersion とよぶ。

平均と φ の間の関係を推定するために、正規化後のデータを使って平均と φ をまず求める。このときの φ は正規化後のデータから平均と分散を \(V = \mu + \mu^{2}\phi\) に代入して計算したものとなる。この φ は正確とは限らないが、大まかな傾向を把握できる。この初期の平均と φ の関係を確認し、φ が平均に依存しないと仮定したときの φ の値を計算する。このときの φ は各遺伝子の平均発現量に依存しないので、すべての遺伝子の平均と分散の関係が、この 1 つの φ でモデル化できる。この φ がすべての遺伝子に共通することから、common dispersion と呼ばれている。次に、common dispersion をベースとして考え、平均と φ の間の関係を局所回帰を行う。局所回帰によって導かれる φ を trended dispersion とよぶ。最後に、最初に計算した個々の遺伝子の粗い φ に対して、treneded dispersion の重みをかけて補正すると、各遺伝子それぞれに対して φ が求まる。このときの φ を gene-wised dispersion とよぶ。各遺伝子の分散を推定するときに、この gene-wised dispersion が使われる。

検定

平均と分散を推定できると、次に一般化線形モデルを利用して平均に差があるかどうかを検定する。例えば、ある遺伝子におけるストレス処理群と対照群の比較について考えてみる。対照群での発現量の期待値を β₁ とおき、ストレス処理群での発現量の期待値と対照群の差を β₂ とおく。このとき、ストレス処理群と対照群における平均発現量に差がなければ、β₂ = 0 であり、そうでなければ β₂ ≠ 0 である。つまり、平均発現量の比較は、β₂ = 0 かどうかの検定に帰着できる。

β₂ = 0 かどうかを検定する方法について考える。まず、ストレス処理群と対照群の期待値に有意差があると仮定したときに、パラメーター β₁ および β₂ で期待値を表す線形モデルについて考える。対照群の期待値は簡単で、E[control] = β₁ である。次に、ストレス処理群では、β₂ ≠ 0 であることから、E[stress] = β₁ + β₂ である。こうして作られるモデルは、想定されるすべてのパラメーターを使っているため、full model とよぶ。なお、これを行列式で書き表すと、次のようになる。

次に、ストレス処理群と対照群の期待値に差がないと仮定したときのモデルを考えてみる。ストレス処理群と対照群の期待値に差がないので、β₂ = 0 であるから、E[control] = β₁ かつ E[stress] = β₁ + 0×β₂ となる。このモデルでは実質的に β₂ を使っていないので、full model に対して reduced model とよぶ。reduced model を行列式で表すと次のようになる。

ここまで full model と reduced model を作った。次に、この 2 つのモデルのうち最も優れているのはどれなのかについて考えていく。こうした線形モデルではどのモデルが優れているのかを比較するために AIC などが使われるが、一般か線形モデルでは尤度比検定あるいは Wald 検定がよく使われている。仮に、尤度比検定あるいは Wald 検定により、full model の方が優れているとわかったならば、これは β₂ ≠ 0 であることを意味するので、この遺伝子の発現量は処理群と対照群の間で有意差がみられると言える。逆に reduced model の方が優れているとわかったならば、β₂ = 0 を意味するので、処理群と対照群の間で、その遺伝子の発現量に有意差が見られないと言える。

edgeR や DESeq2 パッケージでは、β を調整するための行列（0 または 1 からなる行列）のことをデザイン行列という。デザイン行列を定義すると、そのデザイン行列に基づいて full model が構築される。reduced model を構築するとき、デザイン行列の何列目の要素を 0 にするのかを指定して構築する。最後に検定の時に尤度比検定（edgeR, DESeq2）あるいは Wald 検定（DESeq2）を指定する。

多重検定補正

1 つの遺伝子に対して Type I エラーが 5% のもとで検定を行ったとき、その検定結果に 0.05% の Type I エラーが含まれる。n 個の遺伝子に対して、同様な検定を行うと、その検定結果に (1 - (1 - 0.5)ⁿ )% の Type I エラーが含まれるようになる。n = 200 として実際に計算してみるとほぼ 100% となる。つまり、200 個以上の遺伝子にたいして検定を行ったのであれば、その結果はほぼ 100% エラーが含まれていることになる。このような検定結果はほとんど信用できないので、信頼できるように補正を加える必要がある。補正には 2 通りの方法があり、1 つは Type I エラーの閾値を補正する方法であり、もう 1 つは false discovery rate (FDR) を使って補正する方法である。RNA-Seq のデータ解析では、FDR を行う場合がほとんどである。

多重検定補正を説明するために、次のような例を考える。小石とダイヤモンドが混ざっている状態で、両者を分けることを考える。ただし、分類者は小石かダイヤモンドを直接に見ることができない。そこで、例えば、小石とダイヤモンドの密度を測定し、密度順に並べる。次に、ある定数を閾値として設けて、閾値左側を小石、右側をダイヤモンドとして分類することができる。このとき、小石とダイヤモンドの密度の分布の形によって、きれいに分類できない場合がある。例えば、閾値の左にダイヤモンドが混ざっていたり、閾値の右側に小石が混ざっていたりする。

Type I エラーとは、すべての小石のうち、間違って小石でないと判定した石の割合を表す（帰無仮説が正しいのに棄却してしまったエラー率）。上の図で言えば、Type I エラーは、灰色の分布の面積に占めるオレンジ色の部分の面積の割合である。この Type I エラーに着目すると、これは小石をなるべく小石と判定できるように閾値を設けていると言える。閾値を 0.05 にすれば、閾値右側に存在する小石（ダイヤモンドに間違えられる小石）は 5% となる。閾値を 0.01 にすれば、ダイヤモンドに間違えられる小石がわずか 1% となる。

一方で、FDR は、ダイヤモンドと分類されたものの中に混ざっている小石の割合を表す。上の図で言えば、FDR は、閾値右側にある青色の分布の面積に占めるオレンジ色の部分の面積の割合である。閾値（FDR）を 0.1 にすれば、ダイヤモンドとして分類されたものの中に実は 10% の小石が混ざっていることになる。また、FDR を 0.05 にすれば、ダイヤモンドとして分類されたものの中に 5% の小石が混ざっていることになる。

これらのことを発現変動遺伝子の検定にあてはめて考えると、FDR = 0.10 を閾値として用いた時に、検定で得られた発現変動遺伝子の中に高々 10% の遺伝子は実は発現変動遺伝子ではないことになる。RNA-Seq を利用した発現変動遺伝子の検出では、実際に FDR = 0.05 または FDR = 0.10 のように緩い閾値を用いることが多い。緩い閾値を用いることで候補遺伝子を多く検出し、それらの遺伝子の中からストレスに応答する遺伝子をさらに絞っていくのが目的であるから。また、RNA-Seq 実験の場合、候補遺伝子を見つけてから、さらに定量 PCR 実験などで再検証を行うので、初めから候補遺伝子を厳しく絞る必要はない。