1. サポート
2. Help - 解析理論とケーススタディ
3. RNA-Seqデータの解析手法は、これまで偏りすぎていたのではないか？｜DEG以外の見方の可能性

HuR target 200遺伝子の弱い変化から考える、オミクスデータの使い方

前回のケーススタディ「DEG解析で有意とならなかった遺伝子にも、実験で機能が検証されたものがある」では、HuRが脂肪細胞分化を抑制する仕組みを調べた研究を題材に、DEG解析で有意とならなかった遺伝子にも、実験で機能が検証されたものがあることを紹介しました。

このHuR研究では、RNA-Seqだけで候補遺伝子を探したわけではありません。まずRIP-SeqによってHuRに結合するRNAを調べ、HuR immunoprecipitationへのenrichmentが高い200遺伝子をHuR targetとして選びました。そのうえでRNA-Seqを使い、HuR knockoutによってこれらの遺伝子の発現がどのように変化するかを調べています。著者が注目して実験で検証したINSIG1は、通常のDEG解析であれば検証実験の候補に上がらなかったのではないかと考えられます。

今回は、この研究で行われた解析をもう少し深く掘り下げてみます。

200遺伝子を一つずつ評価したわけではない

一般的なRNA-SeqのDEG解析では、各遺伝子についてFold Change、P値、FDRなどを計算し、変化の大きい遺伝子や統計的に有意な遺伝子を探します。しかし、原著論文でHuR target候補200遺伝子について行われた解析は少し違います。

著者らは、それぞれの遺伝子が有意に変化したかどうかではなく、HuR target候補200遺伝子と、それ以外の遺伝子について、HuR knockout/controlのFold Changeの累積分布を比較しました。そしてKolmogorov–Smirnov testを使い、HuR target群が全体としてknockoutで低下する方向に偏っていることを示しました。

つまり、ここで問われているのは、

「どの遺伝子が有意に変化したか」

ではなく、

「HuRに結合する遺伝子群は、他の遺伝子群とは異なる動きをしているか」

です。

Fold Changeの大きい遺伝子が結果を作っているだけではないか

この解析を見ると、一つもっともな疑問が浮かびます。

HuR target候補200遺伝子のうち、ほとんどはWTとKOで明確な発現差はみられず、ごく一部の遺伝子だけで明確な発現差がありました。このことから、HuR target群の累積分布の違いも、実際にはFold Changeの大きい少数の遺伝子によって作られているだけではないか、という疑いです。

これを調べるため、今回のケーススタディでは、HuR target候補のうちFold Changeの大きい遺伝子を上位から順に除き、そのたびにtarget群とnon-target群の分布を比較しました。

なお、以下は原著論文の解析を完全に再現したものではなく、今回利用できたRNA-SeqデータとHuR targetリストを用いた近似的な再解析です。原著論文では200遺伝子がHuR targetとして扱われていますが、今回の再解析ではFold Changeを計算できた199遺伝子を対象としました。

比較対象としたnon-target群では、KO/WTの平均値は約1.090、中央値は約0.995でした。一方、HuR target群では平均値が約0.996、中央値が約0.987でした。non-target群の平均値は一部の大きなFold Changeの影響を受けていますが、中央値で見ると両群とも1付近にあります。

Fold Changeの大きい遺伝子を除外すると、target群のKO/WT中央値は0.9872から0.9910へと、予想どおり1に近づいていきました。つまり、残ったtarget群では、個々の遺伝子の発現変化はさらに小さくなっています。

一方、target群とnon-target群の分布を比較したKS testのP値は0.076から1.1×10⁻⁴へと小さくなりました。KS testは平均値だけではなく分布全体の違いを検出するため、Fold Changeの大きい遺伝子を除いてtarget群のばらつきが小さくなったこと自体が、P値の低下に影響している可能性があります。そのため、P値が小さくなったことをそのまま「生物学的な差が強くなった」と解釈することはできません。

ただし、大きなFold Changeを示す少数の遺伝子を除いても、target群にはKO/WTが1より小さい遺伝子がやや多く、non-target群との分布の違いも残りました。少数の大きなFold Changeだけが原著論文で観察された傾向を作っていた、という説明では十分ではなさそうです。

「変化が小さい」ことは、「情報がない」ではない？

通常のDEG解析では、変化量が大きい遺伝子やP値の小さい遺伝子から見るのは合理的です。しかし今回の再解析の結果は、それだけでは捉えられない情報が存在する可能性を示しています。

HuR target群では、個々の遺伝子のFold Changeは小さく、しかもFold Changeの大きい遺伝子を除外すると、残った遺伝子群の中央値はさらに1に近づきます。それでも、遺伝子群全体としてみるとnon-target群との分布差が残ります。

前の記事でも触れたとおり、遺伝子発現から最終的な生物学的機能に至るまでには、RNA abundanceだけでなく、translation efficiency、protein stabilityなど複数の段階があります。

RNA量では10%（1.1倍）程度の変化しかなくても、translation efficiencyやprotein stabilityでも同じ方向の小さな変化が乗数的に積み重なれば、最終的なタンパク質量や細胞機能への影響はより大きくなる可能性があります。もちろん逆に、RNA量の変化が別の段階で補償されることもあります。

そのため、RNA-Seqだけで候補遺伝子を順位付けしたときに目立たない遺伝子でも、別の情報、このケースではRNAへの結合を組み合わせることで意味が見えてくることがあるのではないでしょうか。

200個の中には、ほかにも重要な候補があるかもしれない

HuR target候補200遺伝子の中で、INSIG1はRNA-Seq上で特に大きな変化を示した遺伝子ではありませんでした。それでも原著論文では、HuRとの直接的な結合、既知の脂肪細胞分化との関係、その後のqPCR、protein、mRNA stability、functional assayなどを組み合わせることで、その役割が検証されました。

そうであれば、200遺伝子の中にはINSIG1以外にも、RNA-Seq単独では目立たないものの、別の実験情報を組み合わせれば重要性が見えてくる候補が存在する可能性があります。

今回のケーススタディでは、RNA-Seqデータに対して、原著論文で行われたFPKMベースのt検定とは別に、FASTQファイルから再計算したGene Countsを用いて、t検定、DESeq2、edgeRによる解析を行いました。DESeq2とedgeRでは、低Counts領域の処理も行っています。

この論文で採用された解析手法は、あり得た解析の道筋の一つにすぎません。別の研究者であれば、異なる手法を選んでいた可能性もあります。そこで、複数の解析手法を適用し、どの解析の道筋を通ったとしても有意候補として拾われやすい遺伝子かどうかを見るために、P<0.05となった手法数を比較しました。

こうして分けた遺伝子群について、過去の論文から、脂肪細胞の分化に直接結びつく機能検証が行われたことのある遺伝子をAIを使って抽出しました。結果を下の表にまとめます。

※ Student's t-test（FPKM）、Student's t-test（counts）、DESeq2、edgeRの4手法について、P<0.05となった手法数を集計した。

興味深いのは、4種類の解析すべてでP<0.05となった、解析方法を変えても有意候補として拾われやすい遺伝子だけに、機能検証例が集中しているわけではないことです。3手法、2手法、1手法、さらにはどの手法でもP<0.05にならなかった遺伝子の中にも、脂肪細胞分化や脂質代謝に関して実験的に機能が示されたものがあります。

これは、少なくとも今回調べた200遺伝子では、「統計学的に有意になりやすい遺伝子ほど、生物学的に重要である」という単純な関係では説明できないことを示しているのではないでしょうか。

ChIP-Atlasの解析も、別の見方ができるかもしれない

以前、ChIP-Atlasを使って転写因子の結合候補とRNA-Seqの発現変化を比較した解析では、結合候補遺伝子の発現変化には上昇と低下が混在し、散布図や個々のFold Changeを見るだけでは明確な傾向を捉えにくいという結果になりました。

しかし、今回のケーススタディで取り上げたHuR研究を見ると、別の問い方ができます。転写因子の結合候補群と非結合遺伝子群について、発現変化の累積分布を比較したらどうなるでしょうか。

個々の候補遺伝子では変化が小さくても、結合候補群全体として何らかの方向性が存在する可能性があります。CDFで可視化し、Kolmogorov–Smirnov testなどで分布を比較すれば、散布図や個々のFold Changeだけでは見えなかった傾向が見つかるかもしれません。

これは単に新しい統計手法を使うという話ではありません。同じRNA-Seqデータに、別の問いを投げているのです。

これまでの解析手法は間違っていたのか

この記事で、これまでのDEG解析を否定したいのではありません。すべての遺伝子を実験で検証できない以上、P値やFold Changeを使って候補を絞り、変化の大きいものや再現性の高いものから調べるのは合理的です。DESeq2やedgeRも、それぞれ明確に定義された統計的問題を解くために開発されています。

その背景には、「重要な遺伝子は上位候補の中により高い確率で含まれている」、あるいは「P値やFold Changeで選んだ遺伝子群には、生物学的に意味のある遺伝子が濃縮されている」という期待があります。この考え方自体はきわめて自然であり、妥当です。

ただし、それらが広く使われてきたことと、RNA-Seqを研究に利用する方法が十分に確立されたことは別の問題です。現在「標準的」と呼ばれる方法も、もとは誰かが考え、試し、論文として発表した一つのアイデアでした。それが有用だったからこそ広く使われ、ソフトウェアとして実装され、多くの派生手法が生まれてきました。

一方で、その成功によって、「RNA-Seqを解析するとはDEGを求めること」「P値やFDRで遺伝子を絞ることが候補探索の基本である」という考え方が強く定着してきた面もあるのではないでしょうか。今回のケーススタディで取り上げたHuR研究では、RIP-Seqによって生物学的に意味のある遺伝子群を先に定義し、RNA-Seqではその集合全体の挙動を調べています。一般的なDEG解析とはかなり異なる発想ですが、生物学的には自然な問いです。

我々はまだ、オミクスデータの解析方法を知らない

RNA-Seqには、正規化、DEG解析、PCA、クラスタリング、GSEA、pathway analysisなど、多くの確立された解析手法があります。しかし、それらは「オミクスデータから研究上重要な情報をどう取り出すべきか」という問題そのものへの答えではありません。

個々の遺伝子を見るのか、遺伝子群を見るのか。RNA量だけを見るのか、translationやproteinなど別の段階と組み合わせるのか。変化の大きい遺伝子を見るのか、別の実験で定義された遺伝子群に共通する弱い変化を見るのか。どのような問いを立てるかによって、同じデータから見えるものは変わります。

これまでオミクス解析では、統計モデル、正規化法、DEG手法などについて非常に多くの議論が行われてきました。一方で、「この生物学的な問いなら、この遺伝子群をこう見たらどうか」といった発想については、まだ議論の余地が大きく残っています。日常的に仮説を立て、実験し、予想外の結果を解釈しているウェットの研究者の視点を、これまで以上に解析方法の議論へ取り入れることが重要ではないでしょうか。

オミクスデータを研究にどう使うべきかという本質的な議論は、まだほとんど始まってすらいません。 一般的な正解がまだないという認識に立ち、ウェットとドライの研究者が、既存の解析フローにとらわれずにデータの見方そのものを議論していく必要があるのだと思います。