この記事は、Case Study No.421 「事例（1）｜低分散・少サンプルのin vitroデータ」の続きです。
前回のような少サンプル・低分散のin vitroデータとは異なり、 今回は、GSE121212のnon-lesional皮膚20サンプルとlesional皮膚20サンプル、 合計40サンプルからなるbiopsy由来のRNA-Seqデータを用います。 各群20サンプルあることで、遺伝子ごとの群内分散は前回より安定して評価しやすくなり、 偶然極端な低分散に見積もられた遺伝子がt検定で有意になってしまうリスクは小さくなります。
一方で、biopsy由来のデータでは、個体差、組織状態の違い、総リード数の違い、 低発現側の測定限界の違いなどが解析結果に反映されやすくなります。 そのため今回は、サンプル数が増えたことで分散推定が安定しやすくなる一方で、 生体サンプル特有のばらつきやダイナミックレンジの違いが残るデータとして、 edgeR、DESeq2、t検定の結果を比較します。
biopsyデータでは、サンプルごとのダイナミックレンジが大きく異なる
Fig1: GSE121212の40サンプルにおけるGene Counts分布の比較。 正規化後のヒストグラムでは、右側のシグナル領域はそろう一方で、 左側の低Count領域では、サンプルごとにヒストグラムの伸び方が異なります。 下の表では、40サンプル中36サンプル以上でGene Countsがない遺伝子を除去した後、 各サンプルの正規化後Gene Countsをlog2変換し、その25th percentileの値が高い順にサンプルを並べています。
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: edgeR、DESeq2、t検定でp < 0.05となった遺伝子の比較。 多くの遺伝子は3つの方法で共通して有意と判定されました。 一方で、edgeRまたはDESeq2では有意だがt検定では有意でない遺伝子、 および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 中・右）
3手法で共通して有意な遺伝子群は、群間差として妥当に見える
Fig3A: edgeR、DESeq2、t検定の3つすべてでp < 0.05となった遺伝子群のヒートマップ。 左はintersection全体、中央はそのうちfold change 1.4倍より大きい遺伝子群、 右はfold change 1.4倍以下の遺伝子群を示します。 fold change 1.4倍以下の遺伝子群でも、non-lesionalとlesionalの間に一貫した差が見られます。
3つの方法で共通して有意と判定された遺伝子群は、 ヒートマップ上でもnon-lesionalとlesionalをよく分けています。 この結果は、3手法のintersectionが、比較的安定したDEG候補を抽出していることを示しています。
興味深いのは、fold change 1.4倍以下の遺伝子群でも、 群間差として解釈しやすいパターンが残っていることです。（Fig3A 右） 前回の少サンプル・低分散データでは、 p値だけで抽出すると変化量の小さい遺伝子が多く含まれ、 fold change条件で絞り込むことに実用的な意味がありました。
しかし、今回のようにサンプル数がある程度あるbiopsyデータでは、 変動幅が小さくても、多数のサンプルで一貫して観測される差があります。 とくに小さな発現変化をターゲットにする研究では、 fold change条件で機械的に除外せず、p値だけで抽出した遺伝子群を確認する価値があります。
方法特異的な遺伝子群は、発現領域ごとに分けて確認する必要がある
Fig3B: 方法特異的に有意と判定された遺伝子群のヒートマップ。 左はedgeRまたはDESeq2でのみ有意となった低Count領域の遺伝子群、 中央はedgeRまたはDESeq2でのみ有意となったシグナル領域の遺伝子群、 右はt検定のみで有意となった遺伝子群です。 下部のグラデーションは、各サンプルの正規化後Gene Countsにおける25th percentileを示しています。 色が濃いほど25th percentileが高く、ダイナミックレンジが狭いサンプルであることを表します。
Fig3Bでは、3手法共通ではなく、特定の方法でのみ有意となった遺伝子群を確認しています。 ここでは、edgeRまたはDESeq2で有意となり、t検定では有意とならなかった遺伝子を、 non-lesionalとlesionalの両方で平均Gene Countsが20を上回るシグナル領域の遺伝子（約1000個）と、 それ以外の低Count領域の遺伝子（約2000個）に分けて詳しく見ていきます。
Fig3B左の低Count領域では、一見するとlesionalで発現が低下しているように見える遺伝子群があります。 しかし、ヒートマップ下部の25th percentileを見ると、 そのパターンはダイナミックレンジの狭いサンプルの偏りと対応しているように見えます。 この場合、ヒートマップ上の群間差は、生物学的な発現変化ではなく、 サンプルごとのダイナミックレンジの違いを反映している可能性があります。
ここで注意すべきなのは、 ダイナミックレンジの狭いサンプルがどちらの群に偏るかは、 必ずしも生物学的条件によって決まるわけではないという点です。 低Count領域では、測定レンジの偏りが群間差のように見えることがあります。 そのため、edgeRやDESeq2で有意と判定された低Count領域の遺伝子は、 p値だけに頼って判断するのは危険です。
一方、シグナル領域を見ると、edgeRまたはDESeq2でのみ有意となった遺伝子群は対角線の上側、 t検定のみで有意となった遺伝子群は対角線の下側で、 いずれも判定境界付近に並んでいました。 これらは、強い発現差を示す遺伝子群というより、 正規化方法や変化量推定のわずかな違いによって、 p値の境界を越えたり越えなかったりした遺伝子群に見えます。
以上をまとめると、edgeRやDESeq2のみで有意と判定される遺伝子は数が多いものの、その多くは低Count領域にありました。 これらの遺伝子については、p値だけに頼らず、正規化後Gene Countsのヒストグラムや線グラフを使って可視化しながら慎重に確認する必要があります。 一方、シグナル領域では、手法の違いによって抽出される遺伝子の差分の多くが判定境界付近に位置しており、 強い発現差を示す別個の遺伝子群というより、判定境界付近で有意・非有意が分かれた遺伝子群と考えられます。
低Count領域では、測定値の有無に基づくフィルタリングが有効
Fig3C: Fig3B左の低Count領域に対して、ON in lesionalとOFF in lesionalの遺伝子群を重ねて表示した図。ON/OFF型遺伝子は、OFF側では40%以上のサンプルで測定値がなく、ON側では60%以上のサンプルで測定値がある遺伝子として抽出しました。
低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値の閾値付近に位置しており、 正規化、分散推定、p値計算、変化量推定の違いによって、 そのため、3手法共通の遺伝子群を優先し、 方法特異的な遺伝子群は、必要に応じて候補を広げる場合に加えるとよいでしょう。
第四に、低Count領域では、連続的な発現量の差としてp値を計算するよりも、 測定値の有無を基準にしたフィルタリングが有効です。 ON/OFF型の変化を見たい場合には、 一方の群で測定され、もう一方の群で測定されないという条件を直接使う方が、 ダイナミックレンジの差を反映しかねないp値よりも、解釈しやすいことがあります。
標準的な解析プロトコルを考える
ここまでの解析では、3つの方法そのものの違いを比較するために、 edgeRとDESeq2には全遺伝子のGene Countsを入力し、 t検定ではlog2変換とglobal normalizationのみを行った正規化後Gene Countsを用いました。 したがって、ここまでの比較は、 edgeRやDESeq2の推奨ワークフロー全体を評価したものではなく、 低Count領域やダイナミックレンジの違いを含むデータに対して、 各手法がどのような遺伝子を有意と判定しやすいかを確認するためのものでした。
実際の解析では、edgeRやDESeq2でも、低発現遺伝子のフィルタリングを行ったうえで検定に進みます。
edgeRでは、filterByExpr を用いて低発現遺伝子を除外しました。
この関数は、単純な合計Countではなく、ライブラリサイズを考慮したCPMを基準に、
一定数以上のサンプルで十分な発現量を持つ遺伝子を残します。
今回のデータでは各群20サンプルであるため、既定値ではおおよそ17サンプル以上で
基準以上のCPMを持ち、かつ全サンプル合計で一定以上のCountを持つ遺伝子が検定対象になります。
その後、TMM正規化、分散推定、quasi-likelihood GLMによる検定を行います。
一方、DESeq2では事前フィルタリングを行ったうえで、
DESeq と results を実行します。
今回の解析では、DESeq2のvignetteで示されるpre-filteringの例に従い、
少なくとも10カウントを持つサンプル数が最小グループサイズ以上である遺伝子を残しました。
今回のデータでは各群20サンプルであるため、
40サンプル中20サンプル以上で10カウント以上ある遺伝子が検定対象になります。
また、results では independent filtering が既定で行われます。
そして、t検定の前処理はSubio Platform上でデータ分布を確認しながら行いました。 具体的には、正規化後Gene Countsが20未満であれば20に置換し、 欠損値を16で補完し、non-lesionalとlesionalの両群において半数以上のサンプルで Counts値が20未満の遺伝子を解析から除外したうえで、t検定を行いました。
これらの処理は、いずれも低Count領域の不安定性やダイナミックレンジの差に由来する 見かけ上の発現変動に基づくDEG候補の混入を大きく減らすのに貢献します。 その結果、3種類のプロトコルで抽出されたDEG候補は、大きく重なりました。
ただし、手法間の違いも若干ながら見られ、その特徴は前のケーススタディと一致します。 すなわち、t検定は観測されたばらつきが極めて小さい場合に有意判定しやすく、 逆に、ノイズ領域の境界付近でばらつきが急に大きくなってくる領域では、 edgeRやDESeq2は観測されたばらつきだけでなく、 同じくらいのCounts値を持つ他の遺伝子のばらつきを考慮するため、 t検定より有意と判定されやすい傾向が見られました。
Fig4: edgeRおよびDESeq2の標準的な推奨ワークフローに近い解析結果と、 Subio Platform上で前処理を行ったt検定結果との比較。 上段は、p < 0.05 かつ fold change 1.4倍以上で抽出された遺伝子を、 non-lesional群とlesional群の平均Gene Counts上に表示した散布図です。 中段は、edgeR、DESeq2、t検定で抽出された遺伝子の重なりを示しています。 下段は、手法間で判定が分かれた遺伝子の発現パターンを、線グラフで示しています。
Fig4のとおり、低発現遺伝子のフィルタリングや低Count領域の前処理を適切に行うと、 edgeR、DESeq2、t検定の主要なDEG候補は大きく重なります。 したがって、信頼性の高い結果を得るには、 検定方法そのものだけでなく、検定を行う前にデータ分布を確認し、 前処理とフィルタリングを適切に設定することが重要です。
ただし、低Count領域の遺伝子をすべて解析対象から外してしまうと、 一方の群では測定され、もう一方の群では測定されないようなON/OFF型の遺伝子を見落とす可能性があります。 上述のとおり、このような遺伝子群は、p値に頼らず、 測定値の有無に基づいて抽出することができます。 そこで、Subioでは次の手順を推奨します。
Subioが推奨する解析手順
正規化後のGene Countsに、low signal cutoffを設定します。 このcutoffは、あらかじめ「10」「20」「30」のように固定値として決めるものではなく、各データセットの分布を確認したうえで設定します。特に、最もダイナミックレンジが狭いサンプル群において、 シグナル領域として扱える下限付近に設定します。 その下限値より上を、連続的な発現量として比較するシグナル領域、 その下を、測定値の有無に基づいて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を用いる場合には、 ノイズ領域との境界付近で、ばらつきの大きい遺伝子が有意と判定されやすいため、 発現パターンを確認し、必要であればCVなどを指標にフィルタリングして除外します。 一方、t検定を用いる場合は、 p値だけでなくfold change条件を組み合わせることで、 変動幅が小さいにもかかわらず、群内分散が小さいために有意となった遺伝子を除外できます。
いずれの方法を用いる場合でも、 低Count領域のON/OFF型遺伝子は別枠で抽出することで、 p値だけでは扱いにくい候補も安全に拾うことができます。
| 領域 | intersectionを使う場合 | edgeR/DESeq2を使う場合 | t検定を使う場合 |
|---|---|---|---|
| シグナル領域 low signal cutoffより上 |
edgeR、DESeq2、t検定の3手法で共通して有意となる遺伝子をDEG候補とする。 | edgeRまたはDESeq2の標準的なプロトコルに従い、p値を用いてDEG候補を抽出する。ただし、ノイズ領域との境界付近では、ばらつきの大きい遺伝子が有意と判定されやすいため、発現パターンを確認し、必要であればDEG候補から除外する。 | t検定のp値を用いてDEG候補を抽出する。ただし、fold change条件を組み合わせることで、変動幅が小さい遺伝子を除外する。 |
| 低Count領域 low signal cutoffより下 |
p値はそのまま使わず、測定値の有無に基づいてON/OFF型の変化として抽出し直す。 | ||
この記事で使用した解析結果を含むSSAファイルは、
こちら
からダウンロードできます。
SSAファイルを
Subio Platform
に読み込むことで、ヒートマップ、散布図、遺伝子リストなどを実際に操作しながら確認できます。
関連トピック
edgeR・DESeq2・t検定の使い分け｜原理とデータタイプ別の実践ルール
edgeR・DESeq2・t検定を比較した事例（1）｜低分散・少サンプルのin vitroデータ
edgeR・DESeq2・t検定を比較した事例（3）｜paired designで何が変わるのか
edgeR・DESeq2・t検定を比較した事例（4）｜低インプットRNA-Seq（low-input RNA-seq）では何が起きるのか