RNA-Seqの発現差解析では、edgeRやDESeq2のような、 RNA-SeqのGene Countsデータに特化した統計手法がよく使われます。 これらの方法は、Gene Countsデータの特徴を考慮した統計モデルに基づいています。
一方で、実務では、 「RNA-SeqではedgeRを使うべき」 「DESeq2を使うべき」 「RNA-Seqにt検定を使ってはいけない」 といった単純化された説明を見かけることもあります。 この記事では、実際のRNA-Seqデータを用いて、 edgeR、DESeq2、そして前処理後のデータに対するt検定の結果を比較します。
データの特徴：低分散・少サンプルのin vitroデータ
この記事で使用するデータは、 GSE49110 です。 これはbiopsy由来のデータではなく、MDA-MB-231細胞を用いたin vitro実験のRNA-Seqデータです。
比較対象は、ERRαを標的とするsiRNA処理群と、control siRNA処理群です。 各条件は2サンプルずつ測定されているため、サンプル数は各群 n = 2 と非常に少ないです。 一方で、同じ条件内の繰り返しサンプルはかなり安定しており、 群内分散が比較的小さいデータであることが特徴です。
実際の研究データでは、このような少サンプルのin vitro実験は珍しくありません。 そこで、まずこのタイプのデータを使って、 edgeR、DESeq2、t検定でどのような違いが出るのかを確認します。
この記事の目的は、どの方法が優れているかを一律に決めることではありません。 それぞれの方法を適用したときに、 どのような発現パターンが有意と判定され、 どのような遺伝子が有意と判定されにくいのかを確認することです。 これは、解析手法の選択や結果の解釈に役立つ実務的な判断材料になります。
比較した3つの方法
ここでは、同じ siC と siE1 の比較について、次の3つの方法を比較しました。
- edgeRによる発現差解析
- DESeq2による発現差解析
- 前処理、正規化、フィルタリング後のデータに対するt検定
edgeRとDESeq2は、RNA-SeqのCountデータを扱うために設計された統計モデルを用います。 一方、t検定は、Raw Gene Countsに直接適用したわけではありません。 また、測定された全遺伝子にそのまま適用したわけでもありません。 低Count領域の不安定な遺伝子をあらかじめ減らし、 一般的な統計手法を適用しやすい状態に整えたうえで実行しました。
t検定に用いた前処理、正規化、フィルタリングの手順は次のとおりです。
- Log2変換を行う。この時点で0は欠損値になります。
- 70th percentileを用いてglobal normalizationを行う。
- 正規化後、linear scaleで20未満の値を20に置換する。
- 欠損値をlinear scaleで16として補完する。
- siC controlサンプルの平均値に対する比に変換する。
- すべてのグループで少なくとも1サンプル以上Gene Countsが20未満となる遺伝子を除外する。
p値だけで判定すると、手法間の違いは大きく見える
まず、edgeR、DESeq2、t検定のそれぞれで p < 0.05 となった遺伝子を比較しました。
個別サンプル同士の散布図では、低Count領域で測定値が少し違うだけでも、 fold changeとしては大きく見えることがあります。 そのため、対角線から離れた点が広がり、低Count領域のばらつきが大きく見えます。
一方、Fig1の散布図では、siC群2サンプルの平均値と、 siE1群2サンプルの平均値を比較しています。 群平均で比較すると、個別サンプル間のランダムなばらつきはある程度平均化されます。 そのため、低Count領域の点は、個別サンプル同士を比較した場合ほど大きく広がりません。
しかし、群平均の点が低Count領域で対角線付近に集まって見えるからといって、 その領域の測定が安定しているわけではありません。 むしろこの領域では、真の発現差よりもランダムな測定揺らぎの影響が大きくなりやすく、 生物学的に意味のある発現差として解釈しにくい値が多く含まれます。
edgeRやDESeq2は、Gene Countsデータ全体から平均発現量と分散の関係を推定し、 その統計モデルに基づいて有意性を評価します。 しかし、低Count領域を中〜高Count領域から連続する分散トレンドの一部として扱うと、 シグナルよりもノイズの影響が大きい遺伝子まで、有意性評価の対象に入ることがあります。 （Fig1 左・中央）
つまり、edgeRやDESeq2を用いる場合でも、 低Count領域をそのまま解釈対象に含めてよいとは限りません。 実際の解析では、これらの方法と組み合わせて低発現遺伝子のフィルタリングがよく行われます。 このことは、低Count領域が統計モデルに任せて安全に解釈できる領域ではなく、 解析者があらかじめ、どこまでを解析対象に含めるかを判断すべき領域であることを示しています。
Fig1: edgeR、DESeq2、t検定で p < 0.05 となった遺伝子の分布比較。 黒い点は全遺伝子を示し、色付きの点は、それぞれの方法で有意と判定された遺伝子を示します。 3つの方法で有意遺伝子の全体的な分布は似ていますが、 低発現領域や判定境界付近では違いが見られます。
edgeRとDESeq2の違いは、主に判定境界付近に見られる
Fig2を見ると、3つの方法で共通して有意と判定された遺伝子が多い一方で、 方法によって有意・非有意の判定が分かれた遺伝子も見られます。
edgeRとDESeq2の違いは、主に有意・非有意の境界付近に見られました。 （Fig2 左） 明らかに大きな差を持つ遺伝子では、2つの方法の結果はよく一致します。 しかし、p値が閾値付近にある遺伝子では、 どちらの方法を使うかによって、有意と判定されるかどうかが変わります。 今回のデータでは、edgeRの方が低Count領域でも多くの遺伝子を有意と判定する傾向が見られました。
少サンプルのRNA-Seqデータでは、edgeRがDESeq2より多くの候補遺伝子を検出することがあると報告されています。 今回の比較でも、edgeRはDESeq2より多くの遺伝子を有意と判定しており、 その傾向と一致します。 ただし、これはedgeRが常に多くの遺伝子を検出するという意味ではありません。
Fig2: edgeR、DESeq2、t検定で p < 0.05 となった遺伝子の比較。 下段は、3つの方法で検出された遺伝子の重なりを示しています。 ベン図でグレーに重ねた領域の遺伝子を、上段の散布図で黒い点として表示しています。
edgeRとDESeq2のみで検出された遺伝子は、平均値の差がやや大きいことが多い
t検定のみで検出された遺伝子と比べると、 edgeRやDESeq2でのみ検出された遺伝子は、対角線からやや離れた位置にあります。 （Fig2 中央） また、これらの遺伝子群には、群内のばらつきが比較的大きく見える遺伝子も含まれます。 （Fig2B 中央）
これは、edgeRやDESeq2が、各遺伝子の4サンプルだけで観測されたばらつきに依存しているわけではなく、 同じくらいのGene Countsを持つ多くの遺伝子から推定された分散トレンドも利用しているためと考えられます。
今回のデータは全体として群内分散が小さいため、 同じくらいのGene Countsを持つ遺伝子から推定される分散トレンドも小さくなりやすいと考えられます。 その結果、個別の測定値を見るとややばらつきが大きく見える遺伝子でも、 edgeRやDESeq2のモデル上では有意と判定された可能性があります。
t検定のみで検出された遺伝子は、平均値の差が小さいことが多い
一方、t検定のみで検出された遺伝子の多くは、対角線に非常に近い位置にあります。 （Fig2 右） また、これらの遺伝子は群内分散が極めて小さいパターンを示します。 （Fig2B 右）
繰り返し数が2サンプルしかない場合、観測された値は母集団分布の中でも高確率で起こる範囲に入りやすくなります。 その結果、観測された分散が、実際の母集団分散よりもかなり小さく見積もられることがあります。 つまり、群平均の差が小さくても、たまたま群内ばらつきが極端に小さく観測されると、 その遺伝子は有意と判定されることがあります。
p値だけで判定する限り、この結果は、 「RNA-Seqデータにはt検定ではなく、edgeRやDESeq2を用いるべきだ」 という主張を支持しているように見えます。
Fig2B: 3つの方法で共通して検出された遺伝子群、 edgeRとDESeq2のみで検出された遺伝子群、 t検定のみで検出された遺伝子群の発現パターン比較。 下段のベン図でグレーに重ねた領域の遺伝子を、上段の線グラフで黒線として表示しています。
前処理とフィルタリングは、低Count領域に対するリスク管理である
前処理やフィルタリングを行わずにt検定を適用すると、 低Count領域の遺伝子が有意と判定され、DEGリストに混入することがあります。 この領域は、Count値の離散性、欠損値、低シグナル由来の不安定性の影響を受けやすく、 p値だけで発現差を解釈するには注意が必要です。
そこで、正規化とlog2変換のみを行ったデータで有意となった遺伝子と、 Low Signal Cutoff、欠損値補完、フィルタリングを行った後に有意となった遺伝子を比較しました。
Low Signal Cutoff、欠損値補完、フィルタリングを行うと、 低Gene Counts領域でt検定により有意と判定される遺伝子は減少しました。 この結果は、主要なリスクがt検定そのものではなく、 低Count領域の不安定なデータに対して統計検定を直接適用することにあることを示しています。
edgeRやDESeq2は、RNA-Seq Countデータ用の統計モデルを自動的に適用します。 これは非常に便利です。 しかし、統計モデルの内部では、解析者が生物学的な目的に応じて、 どのCount領域までを解釈対象に含めるかを柔軟に調整する余地は限られています。 Fig1の左と中央の散布図を見ると、各群2サンプルという条件に対して、 有意判定はやや積極的に見えます。
edgeRやDESeq2を使う場合でも、低発現遺伝子に対する何らかのフィルタリングと組み合わせる方が適切です。
一方、データを見ながら手作業で前処理とフィルタリングを行う場合には、 Low Signal Cutoffの閾値を解析目的に応じて柔軟に調整できます。 探索的解析では低めの閾値を許容できることがあります。 一方、より厳密なbiomarker候補抽出では、高めの閾値を設定する方が適切なことがあります。 解析目的に応じてこの判断を変えられることには、生物学的な解釈上の意味があります。
Fig3: t検定に対する前処理とフィルタリングの効果。 左は正規化とlog2変換のみを行ったデータにt検定を適用した結果です。 中央は、さらにLow Signal Cutoff、欠損値補完、フィルタリングを行った後の結果です。 右は、前処理とフィルタリングによってt検定の有意遺伝子リストから除外された遺伝子の分布を示しています。 除外された遺伝子の多くは低Gene Counts領域にあり、 前処理とフィルタリングによって、不安定な低Count遺伝子が有意と判定されるリスクを減らせることが分かります。
p値にfold change条件を加えると、DEGリストは実務的に解釈しやすくなる
実際の発現差解析では、p値だけでDEGリストを作るとは限りません。 生物学的に解釈しやすい変化に絞るために、fold change条件を加えることがよくあります。 ここでは、edgeR、DESeq2、そして前処理・正規化・フィルタリング後のt検定について、 p < 0.05 かつ fold change 1.4倍以上という条件で比較しました。
p値だけで判定すると、平均値の差が小さい遺伝子でも有意となることがあります。 とくにt検定では、群内分散が極めて小さい遺伝子は、 fold changeが小さくても有意と判定されることがあります。
fold change条件を加えると、対角線に近い小さな変化の遺伝子が除外され、 有意遺伝子は、より大きな発現差を示す領域に絞られます。
Fig4: edgeR、DESeq2、t検定で p < 0.05 かつ fold change 1.4倍以上となった遺伝子の分布比較。
fold change条件を加えると、手法間の違いは小さくなる
Fig5では、p < 0.05 かつ fold change 1.4倍以上という条件を加えたときに、 3つの方法で得られるDEGリストがどのように変化するかを比較しています。
p値だけで比較した場合、t検定のみで有意と判定された遺伝子は355個ありました。 これらの多くは、群内分散が極めて小さい一方で、平均値の差は小さい遺伝子でした。 また、edgeRやDESeq2のみで検出された遺伝子の中にも、 有意・非有意の判定境界付近にあり、発現変化が小さい遺伝子が多く含まれていました。
fold change 1.4倍以上という条件を加えると、 このような小さな発現変化の遺伝子の多くが除外されます。 その結果、edgeR、DESeq2、t検定で得られるDEGリストの重なりは大きくなり、 p値だけで比較した場合よりも、手法間の違いは小さくなりました。 （Fig5）
RNA-Seqの発現差解析では、p値とfold changeを組み合わせてDEGリストを作ることがよくあります。 したがって、この結果は、実務的なRNA-Seq発現差解析では、 手法間の違いが小さくなることを示しています。
Fig5: p < 0.05 かつ fold change 1.4倍以上で抽出したDEGリストの比較。 fold change条件を加えることで、p値だけで比較したときに見られた手法間の違いが小さくなり、 3つの方法で得られるDEGリストの重なりが大きくなります。
ただし、fold change条件を加えても、手法間の差は完全には消えない
fold change 1.4倍以上という条件を加えると、 t検定のみで有意となっていたわずかな変動幅の遺伝子数は大きく減少します。 一方で、edgeRやDESeq2のみで検出される遺伝子はまだ残ります。 （Fig5）
これらの遺伝子は、群内分散が比較的大きいために、 t検定では有意と判定されなかったと考えられます。 （Fig5B 左・中央） 一方、edgeRやDESeq2は、Gene Countsデータ全体から推定した分散トレンドも用います。 そのため、個別の発現プロファイルを見るとばらつきが大きく見える遺伝子でも、 モデル上では有意と判定されることがあります。
この結果は、fold change条件を加えても、 手法間の違いが完全には消えないことを示しています。 edgeRやDESeq2で得られたDEGリストについても、 個々の遺伝子の発現パターンを確認する必要があります。
Fig5B: p < 0.05 かつ fold change 1.4倍以上で抽出した遺伝子群の発現パターン比較。 fold change条件を加えた後でも、 edgeRやDESeq2のみで検出された遺伝子には、個別サンプルを確認すると群内分散が比較的大きいものが含まれます。
この結果は、低分散・少サンプルのin vitroデータに強く依存している
ただし、今回使用したデータはin vitro実験データであり、 繰り返しサンプル間のばらつきが非常に小さいデータです。 また、各条件の繰り返し数は2サンプルしかなく、 統計的には極めて少ないサンプル数です。
このような、繰り返し数が非常に少ないin vitro実験は、 実際のRNA-Seq研究でも珍しくありません。 低分散・少サンプルのデータでは、 Gene Countsデータ全体から推定される分散トレンドも小さくなりやすくなります。 その結果、t検定では有意にならない遺伝子でも、 edgeRやDESeq2では有意と判定されることがあります。 これは、このタイプのデータでedgeRやDESeq2を使う際に理解しておくべき特徴の一つです。
次の記事では、 biopsy由来の中規模RNA-Seqデータ を用いて、edgeR、DESeq2、前処理後のt検定を比較します。 今回見られた傾向が、どの程度再現されるのかを確認します。
重要なのは、どの手法を信じるかではなく、データを見て判断すること
今回の比較では、edgeR、DESeq2、t検定の結果は大きく矛盾しているわけではありません。 p値だけで比較すると手法間の違いは大きく見えます。 しかし、fold change条件を加えると、 実務的なDEG候補リストの違いはかなり小さくなりました。
一方で、方法特異的な差分に含まれる遺伝子を個別に確認すると、 それぞれ異なる特徴が見えてきます。 t検定のみで検出された遺伝子は、 群内分散が極めて小さい一方で、平均値の差が小さいものが多く含まれていました。 これらの遺伝子は、fold change条件を加えることで大きく減少しました。
逆に、edgeRやDESeq2では有意となり、t検定では有意とならなかった遺伝子には、 個別サンプルを確認すると群内分散が比較的大きく見えるものが含まれていました。 低分散・少サンプルのin vitroデータでedgeRやDESeq2を使う場合には、 このような遺伝子がDEGリストに含まれることを理解し、慎重に解釈する必要があります。
Gene Countsを用いたRNA-Seqの発現差解析では、 edgeRやDESeq2のようにカウントデータ全体の分散構造を考慮できる方法を基本にするのが安全です。 一方で、t検定も、正規化・前処理・フィルタリング・fold change条件・可視化確認を組み合わせれば、 edgeRやDESeq2と大きく矛盾しない実用的な発現差評価方法として利用できます。
統計手法は、候補遺伝子を絞り込むための重要な道具です。 しかし、最終的にどの遺伝子を解釈対象にするかを判断するには、 p値だけでなく、実際のデータを確認する必要があります。
Subio PlatformでDEGリストの違いを可視化する
RNA-Seqデータ解析チュートリアル では、データの取り込み、正規化、フィルタリング、発現差解析の方法だけでなく、 各ステップで何を確認し、どのように判断するかを解説しています。
DEG解析というと、Gene CountsをedgeRやDESeq2に入力し、 ツールに処理を任せるだけだと考える方も多いかもしれません。 しかし、その方法だけでは見落としやすい重要なポイントがあります。
チュートリアルページでは、SSAファイルのダウンロードリンクも用意しています。 SSAファイルはSubio Platformに読み込むことができ、実際のデータを自分で操作しながら確認できます。 この記事で示したようなDEGリストの違いも、ぜひデータを見ながら確認してみてください。
関連トピック
edgeR・DESeq2・t検定の使い分け｜原理とデータタイプ別の実践ルール
edgeR・DESeq2・t検定を比較した事例（2）｜biopsy由来の中規模データ
edgeR・DESeq2・t検定を比較した事例（3）｜paired designで何が変わるのか
edgeR・DESeq2・t検定を比較した事例（4）｜低インプットRNA-Seq（low-input RNA-seq）では何が起きるのか