現在のRNA-Seq解析は、下流の変化を見すぎている?
RNA-Seq解析では、DESeq2やedgeRを用いた差次的発現解析(DEG解析)と、その後のGO解析やPathway解析などのエンリッチメント解析が広く使われています。これらは、条件間で大きく発現が変化した遺伝子や、その遺伝子群に共通する生物学的機能を調べるうえで非常に有効です。しかし最近では、大きく発現が変化した遺伝子だけを見ていても、発現変化を引き起こした原因に近い遺伝子までたどり着きにくいという問題意識も強くなってきています。
2026年に Briefings in Bioinformatics に掲載されたレビュー Multimodal bioinformatic analyses of genome-scale expression beyond gene-centric differential expression では、従来のgene-centricな差次的発現解析を超えて、gene co-expression network、gene regulatory network、multi-omics統合、machine learningなどを利用し、遺伝子発現をより広い制御関係の中で捉える解析の発展が整理されています。レビューでは、単一遺伝子のDE解析だけでなく、ネットワークや複数のomics modalityを組み合わせたsystems-levelの解析へ広がっていることが示されています。
少し雑な言い方かもしれませんが、DEGだけでなく、DEGとして選ばれない遺伝子(non-DEG)や、予想とは逆方向に発現している遺伝子にも注目する流れが広がっていると言えるでしょう。
原因に近い遺伝子を探すのは、そう簡単ではない
注目している生物学的現象の原因に近い遺伝子を探そうとしても、現実にはそう簡単ではありません。
以前のケーススタディでは、ChIP-Atlasを利用し、RNA-Seqで変化した遺伝子群の近くに結合ピークを持つ転写因子から、上流制御因子候補を探索しました。ところが、候補となった転写因子が予想どおりの発現変動を示すことはむしろまれで、候補遺伝子全体が増加または減少のどちらか一方向に偏ることさえありませんでした。これは、原因に近い遺伝子を推定するアルゴリズムの開発が必要である一方で、現実には非常に難しいチャレンジだということを示しています。
そこで、より現実的なアプローチとして、RNA-Seqとは異なる種類の実験データを組み合わせるmulti-omicsの統合解析があります。前々回と前回のケーススタディでは、RNA-SeqとRIP-seqを組み合わせた研究を取り上げました。RIP-seqでは、RNA結合タンパク質(RNA-binding protein: RBP)に結合しているRNAを免疫沈降し、そのRNAをシーケンスすることで、特定のRBPと関連するRNA群を調べることができます。
ここで興味深かったのは、RIP-seqで定義されたRNA群について、個々の遺伝子を見ると通常のDEG解析の閾値を超えないような小さな発現変化であっても、遺伝子群全体として見ると、コントロール集団との分布の違いをKS testで有意な偏りとして検出できる可能性が示されたことです。
詳細は、 「DEG解析で有意とならなかった遺伝子にも、実験で機能が検証されたものがある」 「RNA-Seqデータの解析手法は、これまで偏りすぎていたのではないか?|DEG以外の見方の可能性」 をご覧ください。
今回はmulti-omicsとネットワーク解析を組み合わせる
今回は、少し違うアプローチを試してみたいと思います。
前々回・前回で取り上げた GSE124280 と同様に、RNA-SeqとRIP-seqの両方が取得されている別の研究データ、 GSE86590 を使います。
ただし今回は、
multi-omicsによる実験データと、アルゴリズムによるネットワーク解析を組み合わせる
という方法を試してみます。
使用するネットワーク解析手法は DENetwork です。DENetworkは2025年に Nucleic Acids Research に発表された解析手法で、従来のDE解析では見落とされやすいnon-DE genesにも注目し、遺伝子間ネットワーク上のinformation flowへの影響から候補遺伝子を評価します。論文では、DESeq2などのDE解析は発現量の変化そのものを中心に評価するため、発現変化が小さくてもネットワーク上で重要な役割を持つ遺伝子を見落とす可能性があることが、DENetwork開発の動機の一つとして挙げられています。
ただし、ここでDENetworkを「原因遺伝子を特定する方法」と考えるわけではありません。今回確認したいのは、YBX2のノックアウトという原因側の操作と、DEGとして観察された発現変化という結果の間を埋める候補として、DENetworkがPPIネットワークから選んだ遺伝子を、RIP-seqという別の実験データと照合するとどうなるのかという点です。
使用するデータ:YBX2 knockout brown adipose tissue
GSE86590 では、brown adipose tissue(BAT)におけるRNA結合タンパク質 YBX2 の役割が調べられています。wild typeとYbx2 knockoutマウスについて、室温およびcold exposure後のBAT RNA-Seqが登録されています。さらに、YBX2に結合するRNAを調べるためのYBX2 RIP-seqとIgG controlも取得されています。
原著論文 では、RNA-SeqとRNA immunoprecipitationを統合することで、Pgc1αを含むYBX2のmRNA targetが同定され、YBX2がcold-induced BAT activationにおけるRNA stabilityの調節に関与することが示されています。
つまり、このデータセットには、
- YBX2を欠損させたときのRNA-Seq
- YBX2に結合しているRNAを測定したRIP-seq
という、今回の目的に適した2種類のデータがあります。
DENetworkでYBX2とDown-DEGをつなぐ
今回はsourceを YBX2のみ とし、targetをcold exposure後にYbx2 KOで低下したDown-DEGとしました。つまり、YBX2からDown-DEGまでをPPIネットワーク上で結ぶ経路を探索し、その間に存在するinternal nodeを評価します。
K=3(sourceとtargetを含めて最大3 nodes、つまりinternal nodeを最大1個挟む設定)では該当するpathは得られませんでした。
K=4(internal nodeを最大2個挟む設定)では、
- YBX2 → PCBP1 → KCTD13 → FASN
- YBX2 → PCBP1 → TMEM173 → FASN
という2つのpathが得られました。
さらにK=6(internal nodeを最大4個挟む設定)では、次の5つのpathが残りました。興味深いことに、最終的に残ったpath自体はいずれもsourceとtargetの間にinternal node 2個を入れた4 nodesでした。
- YBX2 → PCBP1 → TMEM173 → FASN
- YBX2 → PCBP1 → KCTD13 → FASN
- YBX2 → PCBP1 → PPP1R9B → FASN
- YBX2 → PCBP1 → GRN → FASN
- YBX2 → PCBP1 → STAU1 → FASN
internal nodeは、PCBP1, TMEM173, KCTD13, PPP1R9B, GRN, STAU1の6遺伝子です。PCBP1はすべてのpathに共通して含まれています。もしこの結果をさらに深掘りするなら、まずPCBP1が今回使用したPPIネットワーク上でhubになっていないかを確認したいところです。PPI上で接続数の多い遺伝子は、多くのpathに含まれやすく、ネットワーク解析の結果にバイアスを与える可能性があるからです。
ここであがってきた遺伝子の発現変化を確認したところ、YBX2ノックアウトでは、Ybx2自身およびFasnでは大きな発現低下(それぞれ0.13倍と0.04倍)が見られました。一方で、internal nodeとして選ばれた遺伝子では、
| Gene | KO / WT Fold Change | p value |
|---|---|---|
| PCBP1 | 0.904 | 0.0051 |
| KCTD13 | 1.323 | 0.0199 |
| TMEM173 | 2.814 | 0.0031 |
| PPP1R9B | 1.316 | 0.0074 |
| GRN | 1.191 | 0.0513 |
| STAU1 | 1.187 | 0.0817 |
となりました。
PCBP1はわずかに減少する一方、KCTD13、PPP1R9B、GRN、STAU1はわずかに増加しています。TMEM173は約2.8倍に増加しています。つまり、YBX2ノックアウトによってDownした遺伝子だけを追っていれば、これら6遺伝子はいずれも候補として拾われません。 DENetworkでは、こうした発現変化の方向や大きさだけでは選ばれない遺伝子が、YBX2とDown-DEGをつなぐinternal nodeとしてあがってきました。
ただし、YBX2 → PCBP1 → TMEM173 → FASNという生物学的経路が証明されたと解釈してはいけません。DENetworkで使われているPPIは既知のタンパク質間相互作用を使ったネットワークであり、このpathそのものが今回のBATで実際に一連の生物学的経路として機能していることを示しているわけではないからです。
DENetworkで選ばれた6遺伝子をRIP-seqで確認する
YBX2 RIP-seqとIgG RIP-seq(コントロール)を比較しました。すると、DENetworkのinternal node 6遺伝子すべてが、YBX2 RIP FPKM >1 かつ YBX2/IgG の比が2倍以上という条件を満たしました。
より厳しい条件では、
- YBX2/IgG の比が2倍以上:6/6
- YBX2/IgG の比が4倍以上:5/6
- YBX2/IgG の比が8倍以上:3/6
となりました。
TMEM173ではIgG側のRIP FPKMが0であったため、YBX2/IgG比そのものは計算できません。一方、YBX2 RIPではFPKM 1.53が検出されているため、今回の閾値による分類では便宜上「8倍以上」の群に含めました。
これは一見すると、DENetworkで選ばれた遺伝子とYBX2 RIP-seqの結果がかなりよく重なっているようにも見えます。しかし、ここで背景を確認する必要があります。DENetworkの解析対象となったPPI/RNA-Seq共通の8,576遺伝子を解析対象全体とし、6つのinternal nodeと、それ以外の遺伝子を比較しました。
FPKM >1を条件とした場合、
| RIP-seq条件 | 背景 | Internal nodes | Fisher片側 p |
|---|---|---|---|
| RIP FPKM >1 & YBX2/IgG比 2倍以上 | 65.3% | 6/6 = 100% | 0.078 |
| RIP FPKM >1 & YBX2/IgG比 4倍以上 | 51.1% | 5/6 = 83.3% | 0.120 |
| RIP FPKM >1 & YBX2/IgG比 8倍以上 | 35.1% | 3/6 = 50.0% | 0.355 |
internal nodeでは確かにRIP enrichmentを示す遺伝子と重複が多く見られましたが、背景でもYBX2 RIP enrichmentを示す遺伝子がかなり多いため、この6遺伝子だけでは統計的に有意なenrichmentとは言えないという結果です。
一方、YBX2 RIP FPKM >5という、より強いRIP signalを要求すると、
| RIP-seq条件 | 背景 | Internal nodes | Fisher片側 p |
|---|---|---|---|
| RIP FPKM >5 & YBX2/IgG比 2倍以上 | 41.7% | 5/6 = 83.3% | 0.049 |
| RIP FPKM >5 & YBX2/IgG比 4倍以上 | 31.6% | 4/6 = 66.7% | 0.084 |
| RIP FPKM >5 & YBX2/IgG比 8倍以上 | 19.4% | 2/6 = 33.3% | 0.331 |
となりました。
RIP FPKM >5かつYBX2/IgG比が2倍以上ではp値が約0.049となりますが、internal nodeはわずか6遺伝子です。また、今回は複数のcutoffを比較しており、このp値について多重比較の補正を行っているわけでもありません。この結果だけを強い統計的証拠として扱うべきではありません。そもそも、RIP-seqでは「FPKMが何以上なら十分なbinding signalとするか」という普遍的なcutoffがあるわけでもありません。そのため今回は、RIP FPKM >1と >5の両方を示し、cutoffによって結果がどの程度変わるかを確認しました。
直感的には、YBX2/IgG比が2倍以上という条件より4倍以上、4倍以上より8倍以上と、条件を厳しくするほど、DENetworkのinternal nodeに強く濃縮されるのではないかと考えたくなりますが、結果は逆でした。例えばRIP FPKM >1の場合、条件を2倍以上から8倍以上へ厳しくすると、背景で条件を満たす遺伝子は65.3%から35.1%へ減少します。しかし、internal nodeでも6/6から3/6へ減少します。背景だけが大きく減ってinternal nodeが残る、という結果にはならなかったため、統計的な差はむしろ弱くなりました。
今回の結果からは、DENetworkが特に「最も強いRIP enrichmentを示すRNA」を選んでいるとは言えません。ただし、これは必ずしも意外ではありません。RIP-seqが見ているのはYBX2とRNAの結合であるのに対し、DENetworkで使っているPPIはタンパク質同士の相互作用だからです。両者は異なる生物学的関係を見ているため、RIP enrichmentの強さとDENetworkの結果が単純に一致するとは限りません。
むしろ重要なのは、異なる原理で得られた2種類の情報が、同じ候補遺伝子でどの程度重なるかを見ることで、機能的に重要な遺伝子を絞り込めるかどうかです。本記事では、その候補遺伝子について生物学的な機能実験を行うことはできないため、本当に機能的に重要かどうかまでは検証できません。それでも、DEG解析だけでは拾いにくい候補を探す方法の一つとして、このような組み合わせを検討する価値はあるでしょう。
アルゴリズムかmulti-omicsか、ではない
RNA-Seq解析で原因に近い遺伝子を探そうとすると、どうしても「どのアルゴリズムを使えばよいか」という話になりがちです。しかし私たちは、RNA-Seq単独のデータ解析だけで原因に近い候補を絞り込もうとするより、マルチオミクスデータと解析手法の組み合わせで探索する方が現実的ではないかと考えています。
もちろん、問題はあります。例えば今回のケーススタディであれば、RIP-seqのcutoffは何が適切か、PPI networkに含まれる実験条件と一致しない相互作用やhub geneへの偏り、解析アルゴリズムのモデルやパラメータへの依存などです。つまり、どのデータもどのアルゴリズムも単独では決定的ではありません。
だからこそ、異なる原理で得られた情報が同じ候補を指すかを見ることに意味があります。そして、その候補をどう解釈し、次に何を検証するべきかについては、異なる専門性を持つ研究者がゼロベースで議論することも重要になるでしょう。
2026年のレビュー が示している「gene-centric differential expressionを超える」という方向も、必ずしもDEG解析を捨てるという意味ではありません。DEG解析を出発点としながら、network、multi-omics、machine learningなど異なる情報を重ねていく方向へRNA-Seq解析が広がっている、と考える方がよいでしょう。