これまでの記事で、RNA-Seqの発現差解析において、 edgeR、DESeq2、t検定で抽出される遺伝子の違いを比較してきました。事例（1） では低分散・少サンプルのin vitroデータ、 事例（2） ではbiopsy由来の中規模データ、そして 事例（3） ではpaired designを考慮した場合を見てきました。
今回は、low-input RNA-Seqのデータを用いて、 入力細胞数が少なくなると、発現差解析の前提となるデータ構造がどのように変わるのかを見ていきます。
使用したデータ
本記事では、GSE130882のGene Countsデータを用いました。 このデータには、異なるinput cell数で調製されたRNA-Seqデータが含まれています。 ここでは、SMART_Nxtのデータを用い、 CD3刺激をコントロールとして、CD3+B7-1 Fc刺激によって発現が変動した遺伝子を、 input cell数 100 cells、1k cells、100k cells で見ていきます。 各条件の繰り返し数は2です。
edgeRおよびDESeq2のp値計算には、raw Gene Countsを用いました。 一方、Subio Platform上での散布図表示、fold change条件の確認、 およびシグナル領域の可視化には、サンプル間のデータ量を補正した正規化後Gene Countsを用いました。
100 cellsでは、条件間比較の前にreplicateの再現性が崩れている
まず、100 cellsと1k cellsで、同一条件内の繰り返しサンプル同士がどの程度一致するかを確認しました（Fig1 左・中）。
Fig1: 100 cellsおよび1k cellsにおけるreplicate比較と条件間平均値比較。 左列はCD3の繰り返しサンプル同士の比較、 中央列はCD3+B7-1 Fcの繰り返しサンプル同士の比較、 右列はCD3とCD3+B7-1 Fcの平均値比較を示します。 散布図は正規化後Gene Counts（Processed Signal）で表示しています。
100 cellsでは、CD3同士、CD3+B7-1 Fc同士のreplicate比較でも、 低〜中シグナル領域に大きなばらつきが見られます。 特に、片方のサンプルでは低い値にとどまり、 もう片方のサンプルでは数百から数千程度の値として見える遺伝子が多数あります。
このような分布は、発現量が連続的に測定されているというより、 少数のRNA分子が逆転写・PCR増幅の初期段階で拾われるかどうかに 強く左右されていることを示していると考えられます。
一方、1k cellsでは、100 cellsに比べてreplicate間の一致が大きく改善しています。 低シグナル領域にはまだばらつきが残りますが、 シグナル領域では連続的な分布が形成され、 比較解析の対象として扱いやすい領域があるように見えます。
100k cellsは高入力条件の参照パターンとして使える
次に、100k cellsのデータを確認しました。 100k cellsでは、同一条件内のreplicate比較がほぼ対角線上に乗っており、 事例（1）の 低分散・少サンプルのデータに近い構造を示しています。
Fig2: 100k cellsにおけるreplicate比較、条件間平均値比較、および3手法のVenn図。 右上の散布図では、CD3とCD3+B7-1 Fcの平均値を比較しています。 Venn図は、100k cellsにおいてedgeR、DESeq2、t検定の各方法でp < 0.05となった遺伝子の重なりを示しています。 中央の6039遺伝子のうち、さらに1.4倍以上の発現変動を示した遺伝子を、 右上の散布図では黒点で示しています。
100k cellsでは、edgeR、DESeq2、t検定の3手法でp < 0.05となった遺伝子の共通部分が大きく、 そこにfold change条件を組み合わせることで、 高入力条件における安定した発現変動候補を定義できます。（Fig2 右）
ここでは、この遺伝子群を「真のDEG」と断定するのではなく、 100 cellsおよび1k cellsの結果を解釈するための高入力リファレンスとして扱います。 通常の実験では、このような高入力条件との比較はできないため、 GSE130882のようなデータはlow-input RNA-Seqの挙動を考えるうえで、非常に参考になります。
100k cellsで得られた参照遺伝子群は、100 cellsと1k cellsでどう見えるか
Fig3では、100k cellsで定義した発現変動候補を黒点として示し、 100 cellsおよび1k cellsの散布図上でどのように見えるかを確認しました。
Fig3: 100k cellsで定義した高入力リファレンス遺伝子群を、 100 cellsおよび1k cellsのCD3/CD3+B7-1 Fc比較上に表示した図。 黒点は、100k cellsで3手法共通かつfold change 1.4倍以上となった遺伝子を示します。 オレンジ色の線は、low-input条件で黒点を比較的高い純度で抽出するために設定した経験的なfold changeおよびシグナル閾値です。
100 cellsでは、100k cellsで発現変動候補として定義された遺伝子群であっても、 発現変動がほとんどない領域から、大きな発現変動を示す領域まで広く散らばっています。 このようなデータから、できるだけ偽陽性が混ざらないようにDEG候補を取り出すならば、 大きめのfold change条件と、高めのシグナル閾値を組み合わせたON/OFF型の抽出条件を用いることが有効と考えられます。（Fig3 左）
一方、1k cellsでは、全体的には対角線に収束する領域が拡大しています。 100 cellsと比べると、1k cellsではCD3+B7-1 Fcによる発現変動を反映した構造が、5000以上の領域では見えるようになっています。（Fig3 右）
しかし、それ以下の大部分の領域では、通常のDEGとして安定して抽出するのは難しいと考えられます。 1k cellsでも、広めのfold change条件と高めのシグナル閾値を組み合わせた ON/OFF型の抽出条件を併用することは有効と考えられます。
p値だけで見ると、100 cellsでは手法ごとに異なる不安定な遺伝子群が抽出される
次に、100 cellsおよび1k cellsについて、 edgeR、DESeq2、t検定でp < 0.05となった遺伝子を、それぞれ黒点で示しました。
Fig4: 100 cellsおよび1k cellsにおいて、edgeR、DESeq2、t検定でp < 0.05となった遺伝子の分布。 黒点は各方法でp < 0.05となった遺伝子を示します。 各パネルのタイトルには、有意と判定された遺伝子数を示しています。 散布図は正規化後Gene Counts（Processed Signal）で表示しています。
Fig5: 100 cellsおよび1k cellsでp < 0.05となった遺伝子と、 100k cellsで定義した高入力リファレンス遺伝子群との重なり。 高入力リファレンス遺伝子群は、100k cellsにおいてedgeR、DESeq2、t検定の3手法で共通してp < 0.05となり、 さらにfold change 1.4倍以上を満たした遺伝子です。
100 cellsでは、edgeRやDESeq2でp < 0.05となった遺伝子が、
連続的なシグナル領域というより、
片方の条件で増幅に乗り、もう片方の条件でほぼ未検出となった領域に多く見られます。（Fig4 上段 左・中）
一方で、t検定では、100 cellsでも比較的広い範囲の遺伝子が有意と判定されます。（Fig4 上段 右）
しかし、いずれの手法も、100k cellsで定義した高入力リファレンス遺伝子群との重なりは3割程度でした（Fig5 上段）。
1k cellsでは、100 cellsに比べてシグナル領域のように見える領域がかなり明確に拡大します。
そのため、一見すると通常のDEG解析が可能に見えます。
実際、p値の計算結果だけを見ると、発現差解析が成立しているようにも見えます。（Fig4 下段）
しかし、Fig3の高入力リファレンスと比較すると、
p値で抽出された遺伝子群と参照遺伝子群との重なりは4〜6割程度でした。（Fig5 下段）
現実の研究データでは、100k cellsに相当する答え合わせ用データは通常ありません。 そのため、このように「それらしく見える」こと自体が解釈上のリスクになります。 つまり、低インプットRNA-Seqのデータでは見た目以上に慎重な解釈が必要です。
以上より、低インプットRNA-Seqデータに統計モデルを当てはめて発現差解析を行うことにはリスクがあります。 解析方針としては、次のような選択が考えられます。
- 広めのfold change条件と高めのシグナル閾値を使ったON/OFF型の条件を組み合わせる。ただし、偽陰性はきわめて多くなります。
- 上記よりも内側まで閾値を緩める。その代わり、偽陽性がかなり多く混ざる可能性を受け入れます。
問題は統計手法ではなく、データ構造である
今回の結果は、edgeR、DESeq2、t検定のいずれについても、 手法そのものが不適切であることを示すものではありません。問題は、特に100 cellsのデータが、通常の発現差検定を適用して解釈できるデータ構造になっていない点にあります。
100 cellsでは、低〜中シグナル領域の値が、 連続的な発現量として安定して測定されているというより、 少数の分子が増幅に乗るかどうかによって大きく変わっているように見えます。したがって、100 cellsのようなultra-low input RNA-Seqでは、 p値に基づく通常のDEG抽出を行っても、 その結果を安定した発現差として解釈することはかなり難しいです。
scRNA-Seq解析にも通じる問題
この結果は、scRNA-Seqのデータ解析で使われる統計手法にも示唆を与えます。 scRNA-Seqでは、1細胞あたりのRNA量がさらに少なく、 観測値は検出・未検出、capture効率、増幅、ライブラリの複雑性、 細胞状態、細胞種構成の影響を強く受けます。
scRNA-Seq専用の統計モデルやワークフローは、 このようなデータ構造を前提として設計されています。 しかし、それらはデータ生成過程の不安定さそのものを消すものではありません。 専用手法を使えば、single-cellデータがbulk RNA-Seqのような 安定した発現量データに変わるわけではありません。 そのため、scRNA-Seqでは、今回のような低インプットのRNA-Seq以上に、 非常に慎重に解釈する必要があります。