注）edgeR全体の説明についてはedgeRを使ってRNA-seqから発言変動遺伝子を見つけるを参考にしてください。
はじめに
RNA-Seqデータ解析におけるedgeRパッケージの分散推定ステップ、estimateCommonDisp、estimateTrendedDisp、そしてestimateTagwiseDispについてですね。これらは発現変動遺伝子（DEG）を正確に検出するために非常に重要な処理です。それぞれの関数について、その目的、背景、推定方法、そして役割を詳しく解説します。
背景：RNA-Seqデータと分散推定の重要性
RNA-Seqのカウントデータは、各遺伝子のリードカウント（その遺伝子からシーケンスされたリードの数）として得られます。このカウントデータには以下のような特徴があります。
- 離散性: カウントは整数値を取ります。
- 平均と分散の関係: 一般的に、平均発現量が高い遺伝子ほど分散も大きくなる傾向があります。
- 過分散 (Overdispersion): ポアソン分布では分散が平均と等しくなりますが、実際のRNA-Seqデータでは、生物学的なばらつき（サンプル間の個体差、実験条件の微妙な違いなど）や技術的なばらつき（ライブラリ調製、シーケンスの偏りなど）により、ポアソン分布で期待されるよりも分散が大きくなる「過分散」がしばしば観察されます。
この過分散を適切にモデル化するために、edgeRでは負の二項分布 (Negative Binomial distribution, NB分布) を用います。NB分布は平均 $\mu$ と分散 $\sigma^2$ の関係を、分散パラメータ $\phi$ (dispersion parameter, edgeRでは通常 $\phi$ と表記され、$\alpha = \phi$ として $\sigma^2 = \mu + \phi \mu^2$ と表されることが多い) を用いてモデル化します。この $\phi$ が大きいほど、平均 $\mu$ に対して分散がより大きくなることを意味します。
正確なDEG解析のためには、この分散パラメータ $\phi$ を各遺伝子に対して正確に推定することが不可欠です。しかし、生物学的反復サンプル数が少ない場合（例えば各群3サンプルなど）、個々の遺伝子のデータだけから安定した $\phi$ の値を推定するのは非常に困難です。推定値が不安定だと、統計検定の信頼性が低下し、偽陽性（実際には変動していない遺伝子を変動していると判定）や偽陰性（実際に変動している遺伝子を見逃す）が増加する可能性があります。
そこでedgeRは、遺伝子間で情報を共有することにより、より安定かつ信頼性の高い分散推定を行うための階層的なアプローチを採用しています。それが、共通分散、トレンド分散、タグワイズ分散の3段階の推定です。
1. d <- edgeR::estimateCommonDisp(d)：共通分散の推定
-
目的:
全ての遺伝子に対して、単一の共通の分散パラメータ $\phi_{common}$ を仮定し、その値を推定します。 -
背景・仮定:
最も単純なモデルとして、データセット中の全遺伝子が同じ分散特性（つまり同じ $\phi$ の値）を共有すると仮定します。実際には遺伝子ごとに分散特性は異なりますが、この共通分散はデータセット全体の大まかな過分散の度合いを捉えるための初期推定値として、また後続の推定ステップの基礎として重要です。 -
推定方法:
edgeRでは、条件付き最尤推定 (Conditional Maximum Likelihood, CML) を改良した手法を用いて $\phi_{common}$ を推定します。CML法は、各遺伝子のグループごとの合計カウントが与えられたという条件のもとで、観測されたカウントデータの配置の尤度を最大化する $\phi$ を求めます。この方法は、グループ間の平均発現量の差の影響を受けにくいという利点があります。
具体的には、各遺伝子についてグループごとの総和を固定した上で、観測されたカウントの分割パターンが生じる確率を計算し、全遺伝子にわたってこの条件付き尤度の積を最大にするような $\phi_{common}$ を探索します。
edgeRは、この計算を効率的かつ頑健に行うためのアルゴリズム（例えば、quantile-adjusted CMLなど）を実装しています。 -
役割:
- データセット全体の平均的な分散レベルを把握します。
- 後述するトレンド分散やタグワイズ分散の推定において、出発点や制約、あるいは事前情報の一部として利用されます。
- 特に生物学的反復数が極めて少ない（例えば各群1サンプルなど、この場合は分散推定自体が非常に困難ですが）場合や、遺伝子数が少ない場合に、最も基本的な分散の推定値となります。
-
DGEListオブジェクトへの影響:
推定された共通分散の値は、DGEListオブジェクトdのd$common.dispersionという要素に格納されます。
2. d <- edgeR::estimateTrendedDisp(d)：トレンド分散の推定
-
目的:
遺伝子の平均発現量（通常はlogCPM: log2 counts per millionで評価）と分散パラメータ $\phi$ の間に見られる系統的な関係性（トレンド）を捉え、発現量に依存した分散パラメータ $\phi_{trended}(logCPM)$ を推定します。 -
背景・仮定:
実際には、遺伝子の平均発現量によって分散の大きさが異なる傾向があります。例えば、非常に低発現の遺伝子ではポアソンノイズ（技術的ノイズ）が相対的に支配的で分散が小さく、一方、高発現の遺伝子では生物学的なばらつきがより顕著になり過分散が大きくなる、といったトレンドが見られることがあります。トレンド分散は、このような発現量依存的な分散のパターンをモデル化しようとします。 -
推定方法:
- 遺伝子ごとの一時的な分散推定: まず、各遺伝子に対して何らかの方法で一時的な分散の推定値を計算します。これは、例えば遺伝子ごとの未調整のサンプル分散や、共通分散を用いた場合の尤度プロファイルから導出されることがあります。
- 発現量との関係の可視化: 次に、これらの遺伝子ごとの分散推定値（またはそれに関連する量）を、対応する平均発現量（logCPM）に対してプロットします。
-
トレンドのフィッティング: このプロットに対して、平滑化手法を用いてトレンドラインを当てはめます。
edgeRでは、遺伝子ごとの対数尤度プロファイル（profile likelihood）を平均発現量に対してプロットし、そのトレンドをノンパラメトリックな方法（例えば、スプラインやLOESSに似たアプローチ）で捉えます。具体的には、各遺伝子のlogCPMにおける、その遺伝子のデータから計算される尤度が最大になるような分散値を求め、それらを平均発現量に対して平滑化することで、平均発現量に応じた分散の期待値（トレンド）を推定します。
このプロセスでは、共通分散から得られた情報をガイドとして利用することもあります。
-
役割:
- 共通分散よりも現実に即した、遺伝子の発現量レベルに応じた分散のベースラインを提供します。
- 個々の遺伝子の分散（タグワイズ分散）を推定する際に、より情報量の多い事前情報として機能します。これにより、各遺伝子の分散をその発現量レベルで期待される値に近づけるように調整し、推定の安定性を向上させます。
-
DGEListオブジェクトへの影響:
推定されたトレンド分散の値は、DGEListオブジェクトdのd$trended.dispersionという要素に、各遺伝子の平均発現量（d$AveLogCPM）に対応する形で格納されます。
3. d <- edgeR::estimateTagwiseDisp(d)：タグワイズ分散（遺伝子ごとの分散）の推定
-
目的:
個々の遺伝子（「タグ」とも呼ばれます）に固有の分散パラメータ $\phi_{tagwise}$ を推定します。これが最終的に統計検定（例：exactTestやGLMベースの検定）で使用される分散の値となります。 -
背景・仮定:
共通分散やトレンド分散は遺伝子間で情報を共有することで安定性を高めますが、それでもなお、個々の遺伝子にはそれぞれ固有の生物学的・技術的な分散特性が存在する可能性があります。タグワイズ分散は、この遺伝子ごとの特異性を捉えつつ、反復数が少ないことによる推定の不安定さを克服するために、経験的ベイズ (Empirical Bayes) の考え方を利用します。 -
推定方法:
経験的ベイズ法とは、集団全体から得られる情報（ここでは共通分散やトレンド分散）を「事前分布」として利用し、個々の遺伝子のデータから得られる情報（「尤度」）と組み合わせることで、より安定した「事後推定値」を得る手法です。
edgeRにおけるタグワイズ分散の推定は、具体的には以下のステップで行われます。- 事前情報の決定: 共通分散 ($\phi_{common}$) とトレンド分散 ($\phi_{trended}$) から、各遺伝子に対する事前分布の中心（期待される分散の値）と、その事前分布のばらつき（信頼度）が設定されます。通常、トレンド分散がより情報量の多い事前情報として用いられます。
- 遺伝子ごとの尤度の計算: 各遺伝子について、その遺伝子の観測カウントデータに基づいて、様々な分散の値に対する尤度が計算されます。
-
事後推定（縮小推定）: 事前分布と尤度を組み合わせることで、各遺伝子の分散パラメータ $\phi_{tagwise}$ を推定します。この際、個々の遺伝子データのみから計算される分散推定値（これは不安定である可能性が高い）は、事前情報（トレンド分散や共通分散）の方向に「縮小 (shrinkage)」されます。
縮小の度合いは、以下の要素によって調整されます。- 個々の遺伝子データの情報量: リードカウントが多い遺伝子や、反復サンプル間で一貫したパターンを示す遺伝子は、その遺伝子自身のデータから得られる情報が重視され、縮小の度合いは小さくなります。
-
事前分布の強さ: 全遺伝子数が多いほど、あるいは共通分散やトレンド分散の推定が安定しているほど、事前分布の信頼性が高まり、個々の遺伝子の分散は事前情報により強く引き寄せられます。
edgeRでは、Cox-Reid profile-adjusted likelihood を最大化するアプローチや、重み付けされた尤度を用いる方法で、各遺伝子のタグワイズ分散を頑健に推定します。
-
役割:
- 各遺伝子の特性を考慮しつつ、反復数が少ない場合でも安定した分散推定値を提供します。
- 分散の推定誤差を小さくすることで、DEG検出における偽陽性を抑制し、検出力を向上させます。
- 過分散が小さい遺伝子はポアソン分布に近い分散値を、過分散が大きい遺伝子はより大きな分散値を持ちつつも、極端に外れた値にならないように調整されます。
-
DGEListオブジェクトへの影響:
推定された各遺伝子のタグワイズ分散の値は、DGEListオブジェクトdのd$tagwise.dispersionという要素にベクトルとして格納されます。このd$tagwise.dispersionが、exactTestやglmFitといった後続の統計検定関数で、各遺伝子の分散パラメータとして直接利用されます。
なぜこの3段階（または2段階）が必要なのか？
- 頑健性と精度: 生物学的反PLING数が少ないRNA-Seqデータでは、個々の遺伝子から直接分散を推定すると非常に不安定になります。共通分散で全体の傾向を掴み、トレンド分散で発現量依存性を考慮し、最後にタグワイズ分散で遺伝子ごとの情報を加味しつつ事前情報で安定化させるという段階的なアプローチにより、頑健性と精度のバランスを取っています。
- 情報の借用: 遺伝子間で情報を「借用」する（borrowing information）という統計的な考え方に基づいています。これにより、個々の遺伝子の情報が少なくても、他の遺伝子から得られた知識を利用して、より信頼性の高い推定が可能になります。
実行順序について
通常、estimateCommonDisp $\rightarrow$ estimateTrendedDisp $\rightarrow$ estimateTagwiseDisp の順で実行されます。estimateCommonDisp は、estimateTrendedDisp や estimateTagwiseDisp のための初期値や参照点として機能します。estimateTrendedDisp は、estimateTagwiseDisp のためのより精緻な事前情報を提供します。
ただし、edgeRのバージョンや使用するワークフローによっては、estimateDisp という単一の関数がこれら複数のステップを内部で実行することもあります。ご提示のコードのように明示的に各ステップを呼び出すことで、解析の流れをより細かく制御し、理解することができます。
これらの分散推定ステップを経ることで、edgeRは信頼性の高い発現変動遺伝子の同定を実現しています。特に反復数が限られている実験デザインにおいては、この精巧な分散推定メカニズムがその性能の鍵となっています。