1. サポート
2. Help - 解析理論とケーススタディ
3. RNA-Seq発現差解析でedgeR・DESeq2・t検定を比較した事例（2）｜biopsy由来・中規模データでの検証

前回のような少サンプル・低分散のin vitroデータとは異なり、 今回は、GSE121212のnon-lesional皮膚20サンプルとlesional皮膚20サンプル、 合計40サンプルからなるbiopsy由来のRNA-Seqデータを用います。 各群20サンプルあることで、遺伝子ごとの群内分散は前回より安定して評価しやすくなり、 偶然極端な低分散に見積もられた遺伝子がt検定で有意になってしまうリスクは小さくなります。

一方で、biopsy由来のデータでは、個体差、組織状態の違い、総リード数の違い、 低発現側の測定限界の違いなどが解析結果に反映されやすくなります。 そのため今回は、サンプル数が増えたことで分散推定が安定しやすくなる一方で、 生体サンプル特有のばらつきやダイナミックレンジの違いが残るデータとして、 edgeR、DESeq2、t検定の結果を比較します。

biopsyデータでは、サンプルごとのダイナミックレンジが大きく異なる

RNA-Seqデータでは、正規化によってすべての領域が均一にそろうわけではありません。 正規化前のGene Countsでは、0付近の左端は比較的そろって見えますが、 右側のシグナル領域の位置は総リード数の違いを反映してサンプル間でずれます。

一方、正規化後のGene Countsでは、右側のシグナル領域の山の位置がそろいます。 これは、サンプル間で比較したい主要な発現領域をそろえるという意味では自然な結果です。 しかし、その結果として、左側の低Count領域では、 サンプルごとにヒストグラムの伸び方の違いが見えるようになります。

ダイナミックレンジが狭いサンプルでは、低発現側の測定値が十分に得られていません。 そのため、正規化アルゴリズムやバッチ補正によって分布をそろえようとしても、 測定されていない低発現側の情報を復元することはできません。 このことは、これまでCase Study No. 403 および No. 413 で確認してきました。

今回の記事では、そのような低発現側の測定限界の違いを含むデータに対して、 edgeR、DESeq2、t検定を適用したときに、どのような結果が得られるのかを確認します。

比較条件

今回は、3つの方法で得られるp値を比較するために、 edgeRとDESeq2には全遺伝子のGene Countsを入力し、 それぞれの統計モデルに基づいてp値を計算しました。 ただし、ここではedgeRやDESeq2の推奨ワークフロー全体を評価するのではなく、 t検定を含めた3つの有意判定方法が、どのような遺伝子を拾うのかを比較することを目的としています。

t検定を行うSubio Platformでも、low signal cutoffや欠損値の補完などの 特別なpreprocessingは行わず、 log2変換とglobal normalizationのみを行ったデータに適用しました。

そのうえで、各方法についてp < 0.05となる遺伝子を抽出しました。 変化量の評価には、non-lesionalに対するlesionalのlog ratioを用いました。

したがって、この解析は標準的な推奨プロトコルを示すものではなく、low signal cutoffや欠損値処理を行わない条件で、3つの有意判定方法の違いを確認するための検証です。

edgeR、DESeq2、t検定の結果は大きく重なる

Fig2を見ると、3つの方法で共通して有意と判定された遺伝子が大部分を占めています。 40サンプルのbiopsyデータでは、主要な発現差については、 edgeR、DESeq2、t検定の3つでかなり共通して検出されています。

この点は、前回記事で扱った少サンプル・低分散のin vitroデータとは大きく異なります。 前回記事では、各群の繰り返しが2サンプルしかなかったため、 有意判定された遺伝子群は、群内の2サンプル間のばらつきの大きさに強く影響されていました。 一方、今回は各群20サンプルが含まれているため、 遺伝子ごとの群内分散や群間差の推定が安定しやすくなり、 主要な発現差については、手法間の違いが小さくなったと考えられます。

一方で、方法依存的に検出される遺伝子も残ります。 とくに、edgeRまたはDESeq2では有意と判定されたものの、 t検定では有意と判定されない遺伝子群が目立ちます。 ただし、この遺伝子群の多くは低Count領域に含まれていました。 edgeRまたはDESeq2でのみ有意となった遺伝子は約3,000個ありましたが、 そのうち約2,000個は低Count領域にあり、 シグナル領域にある遺伝子は約1,000個でした。（Fig3B 左・中）

t検定のみで有意となった遺伝子は956個でした。 そのため、シグナル領域に限って見ると、 edgeRまたはDESeq2でのみ有意となった遺伝子数と、 t検定のみで有意となった遺伝子数の差は大きくありません。（Fig3B 中・右） この結果は、今回のように各群20サンプルが含まれるデータでは、 主要な発現差について、手法間の違いが前回記事ほど大きくならなかったことを示しています。

3手法で共通して有意な遺伝子群は、群間差として妥当に見える

3つの方法で共通して有意と判定された遺伝子群は、 ヒートマップ上でもnon-lesionalとlesionalをよく分けています。 この結果は、3手法のintersectionが、比較的安定したDEG候補を抽出していることを示しています。

興味深いのは、fold change 1.4倍以下の遺伝子群でも、 群間差として解釈しやすいパターンが残っていることです。（Fig3A 右） 前回の少サンプル・低分散データでは、 p値だけで抽出すると変化量の小さい遺伝子が多く含まれ、 fold change条件で絞り込むことに実用的な意味がありました。

しかし、今回のようにサンプル数がある程度あるbiopsyデータでは、 変動幅が小さくても、多数のサンプルで一貫して観測される差があります。 とくに小さな発現変化をターゲットにする研究では、 fold change条件で機械的に除外せず、p値だけで抽出した遺伝子群を確認する価値があります。

方法特異的な遺伝子群は、発現領域ごとに分けて確認する必要がある

Fig3Bでは、3手法共通ではなく、特定の方法でのみ有意となった遺伝子群を確認しています。 ここでは、edgeRまたはDESeq2で有意となり、t検定では有意とならなかった遺伝子を、 non-lesionalとlesionalの両方で平均Gene Countsが20を上回るシグナル領域の遺伝子と、 それ以外の低Count領域の遺伝子に分けて詳しく見ていきます。

Fig3B左の低Count領域では、一見するとlesionalで発現が低下しているように見える遺伝子群があります。 しかし、ヒートマップ下部の25th percentileを見ると、 そのパターンはダイナミックレンジの狭いサンプルの偏りと対応しているように見えます。 この場合、ヒートマップ上の群間差は、生物学的な発現変化ではなく、 サンプルごとの測定レンジの違いや実験要因を反映している可能性があります。

一方で、同じFig3B左の下部には、 non-lesionalでは測定されない、またはほぼOFFに近い遺伝子が、 lesionalでONになるようなパターンも見られます。 このようなON/OFF型の変化は、edgeRやDESeq2のようなcount-based methodが得意とする領域です。

しかし、ここで注意すべきなのは、 ダイナミックレンジの狭いサンプルがどちらの群に偏るかは、 必ずしも生物学的条件によって決まるわけではないという点です。 低Count領域では、測定レンジの偏りが群間差のように見えることがあります。 そのため、edgeRやDESeq2で有意と判定された低Count遺伝子であっても、 p値だけに頼って解釈するのは危険です。

Fig3B中央のシグナル領域でも、注意が必要です。 この群は、両群の平均Gene Countsが20を上回るため、 Fig3B左の低Count領域よりは信号として扱いやすい領域です。 しかし、ヒートマップのパターンを見ると、 下部に示した25th percentileのグラデーションと対応しているようにも見えます。 つまり、シグナル領域に入っていても、 サンプルごとのダイナミックレンジの違いに引きずられたパターンが混ざっている可能性があります。

さらにFig2の散布図と合わせて見ると、 edgeRまたはDESeq2でのみ有意となった遺伝子群は対角線の上側、 t検定のみで有意となった遺伝子群は対角線の下側で、 いずれも判定境界付近に並んでいました。 これらは、強い発現差を示す遺伝子群というより、 正規化方法や変化量推定のわずかな違いによって、 p値の境界を越えたり越えなかったりした遺伝子群に見えます。

Fig3B右のt検定のみで有意となった遺伝子群については、 群間で一貫した大きな発現差を示しているというより、 変動幅が小さく、かつ群内分散も小さいように観測された遺伝子が、 t検定の有意判定を偶然通過した可能性があります。 また、この中にも低Count領域でダイナミックレンジの違いを反映しているように見える遺伝子が含まれています。

したがって、Fig3B中央および右に含まれる遺伝子群は、 主要なDEG候補としてそのまま扱うよりも、 3手法共通の遺伝子群とは分けて、慎重に解釈する必要があります。 今回の解析では、少なくとも優先度の高いDEG候補として扱うより、 追加確認の対象とする方が妥当と考えられます。

低Count領域では、測定値の有無に基づくフィルタリングが有効

低Count領域では、edgeRやDESeq2がON/OFF型の変化を検出することがあります。 しかし、その同じ領域には、ダイナミックレンジの偏りによる見かけの差も混ざります。 そのため、低Count領域の遺伝子をp値だけで判定するのは危険です。

そこでFig3Cでは、edgeRやDESeq2のp値ではなく、 測定値の有無そのものを使ってON/OFF型の遺伝子を抽出しました。 具体的には、OFF側では40%以上のサンプルで測定値がなく、 ON側では60%以上のサンプルで測定値がある遺伝子を抽出しています。この条件により、ヒートマップ上でON/OFF型の変化として確認したい遺伝子群をおおむねカバーできています。

このように、連続的な発現量の差はシグナル領域で評価し、 低Count領域では、測定値の有無に基づいてON/OFF型の変化を抽出することで、 シグナル領域と低Count領域の両方を含めたDEG候補リストを作ることができます。

今回の比較から見えること

今回の結果から、RNA-Seqの発現差解析では、 低Count領域とシグナル領域を分けて考えることが重要であると分かります。

第一に、低Count領域では、有意差判定は難しくなります。 この領域は、正規化後Gene Countsのダイナミックレンジ差を反映しやすく、 群間差のように見えても、実際にはサンプルごとの測定レンジや実験要因の影響である可能性があります。 そのため、低Count領域のp値は慎重に扱う必要があります。

第二に、シグナル領域でedgeR、DESeq2、t検定の3つすべてで有意と判定された遺伝子群は、 かなり妥当なDEG候補に見えます。 今回のようにサンプル数がある程度あるbiopsyデータでは、 変動幅が小さい遺伝子であっても、多数のサンプルで一貫した差が見られることがあります。 そのような遺伝子を対象にする場合、fold change条件を機械的に加えない方がよいこともあります。

第三に、シグナル領域であっても、方法特異的な遺伝子群は慎重に解釈する必要があります。 edgeRのみ、DESeq2のみ、t検定のみで有意となる遺伝子群は、 正規化、分散推定、p値計算、あるいはダイナミックレンジ差の影響を受けている可能性があります。 3手法共通の遺伝子群とは分けて考えるべきです。

第四に、低Count領域では、統計学的手法に頼らなくても、 測定値の有無を基準にしたフィルタリングが強力に使えます。 ON/OFF型の変化を見たい場合には、 一方の群で測定され、もう一方の群で測定されないという条件を直接使う方が、 p値よりも解釈しやすいことがあります。

標準的な解析プロトコルを考える

今回の解析では、3つの方法を比較するために、 あえてlow signal cutoffや欠損値の埋め合わせを行いませんでした。 edgeRやDESeq2の標準的な解析プロトコルでは、 低発現遺伝子のフィルタリング、ライブラリサイズ補正、分散推定、統計モデルに基づく検定が行われます。 しかし、これらの処理だけでは、 サンプルごとのダイナミックレンジの違いによって生じる見かけの発現変動を防げません。

そのため、実際の解析では、edgeRやDESeq2の標準的な手順に加えて、 サンプルごとのデータ分布や、正規化後ヒストグラム左側の伸び方の違いを確認し、 正規化後データに残る低発現側の測定限界の違いを、 そのまま発現差解析に持ち込まないようにする必要があります。

具体的には、正規化後のGene Countsにlow signal cutoffを設定します。 ただし、このcutoffはできるだけ低く設定するのではなく、 最もダイナミックレンジが狭いサンプル群でも、 シグナル領域として扱える下限付近に設定します。 その下限値より上を、連続的な発現量として比較するシグナル領域、 その下を、測定値の有無に基づいてON/OFF型の変化を確認する領域として分けて扱います。

このようにシグナル領域と低Count領域を分けずに、 低発現側の測定限界の違いを含むデータをそのまま解析すると、 サンプルごとのダイナミックレンジの違いが群間差として抽出されることがあります。 このような低発現側の測定限界の違いは、 edgeRやDESeq2の統計モデルにも、標準的な解析プロトコルにも明示的には組み込まれていません。 しかし、現実のRNA-Seqデータには存在し、見かけの発現変動を生む要因になります。 これは単なる理論上の懸念ではなく、 Case Study No. 403 で示したように、 実際のRNA-Seqデータでも観察される問題です。

シグナル領域では、edgeR、DESeq2、t検定の3つの有意判定のintersectionを取る方法が、 堅実なDEG候補抽出として使えます。 3手法で共通して有意となる遺伝子は、 少なくとも今回のデータでは、ヒートマップ上でも群間差として妥当な構造を示していました。

3つの方法をすべて実行してintersectionを取るのが難しい場合には、より簡略な方法として、シグナル領域ではedgeRまたはDESeq2を用いる方法と、t検定を用いる方法が考えられます。

edgeRまたはDESeq2を用いる場合には、 p値で抽出された遺伝子の発現パターンが、 ダイナミックレンジの偏りと対応していないかをヒートマップで確認することをお勧めします。 一方、t検定を用いる場合は、p値で有意判定された遺伝子の中から、 変動幅が小さく、かつ群内分散も小さい遺伝子を除外する処理を加えることで、 t検定がedgeRやDESeq2と比べて弱いとされる点を補うとよいでしょう。

いずれの方法を用いる場合でも、 低Count領域のON/OFF型遺伝子は別枠で抽出する必要があります。

領域	intersectionを使う場合	edgeR/DESeq2を使う場合	t検定を使う場合
シグナル領域 low signal cutoffより上	edgeR、DESeq2、t検定の3手法で共通して有意となる遺伝子をDEG候補とする。	edgeRまたはDESeq2のp値を用いてDEG候補を抽出する。ただし、発現パターンがダイナミックレンジの偏りと対応していないか確認する。	t検定のp値を用いてDEG候補を抽出する。ただし、変動幅が小さく、かつ群内分散も小さい遺伝子は除外する。
低Count領域 low signal cutoffより下	p値はそのまま使わず、測定値の有無に基づいてON/OFF型の変化として抽出し直す。