実験で検証された遺伝子は、RNA-Seqでは目立たないことがある
RNA-Seq解析では、DEG解析を行い、P値やFDRの小さい遺伝子を候補として取り上げることがよくあります。P値やFDRは候補を絞るうえで便利ですが、研究上重要な遺伝子が必ずしも上位に並ぶとは限りません。 研究で知りたいのは、どの遺伝子が最も有意だったかということよりも、どの分子が現象の原因に近く、その作用がどのように表現型につながっているかということです。DEG解析の順位だけでは、その位置づけまでは分かりません。
今回は、GSE124280（RNA-binding proteinであるHuRと脂肪細胞分化の関係を調べた研究）を例に考えてみます。
HuRは脂肪細胞分化を抑える
この研究では、HuRを減少させると脂肪細胞分化が促進し、逆にHuRを増やすと分化が抑えられることが示されています。脂肪組織特異的にHuRを欠損させたマウスでも、脂肪量の増加や代謝異常が観察されています。
ここから、HuRがどの遺伝子を介して脂肪細胞分化を抑えているのか、という問いが出てきます。
RNA-Seqでは、HuR欠損後に起きた変化を見る
著者らは、HuR knockoutマウスとcontrolマウスのBAT、eWAT、iWATについてRNA-Seqを行っています。またGSEAを使って、各組織でどのようなpathwayが変化しているかも調べています。adipogenesis、inflammation、browning、myogenesisなど、組織ごとに異なる変化が見つかっています。
ただ、RNA-SeqやGSEAから分かるのは、HuRを欠損させた後に組織全体で何が起きたかです。HuRが直接どのmRNAに作用しているかを知るには、別の実験が必要です。GSEAで見えているpathwayも、HuRの直接標的というより、HuR欠損後に生じた下流の変化を反映していると考える方が自然です。
HuRが直接結合するRNAをRIP-Seqで調べる
HuRはRNA-binding proteinなので、著者らはRIP-Seqを使ってHuRに結合するRNAを調べています。その結果をもとにHuRへの結合が強い遺伝子をrankingし、上位200遺伝子をHuR targetとして扱っています。
この200遺伝子は、RNA-SeqでP値が小さかった遺伝子を集めたものではありません。HuRがどのRNAに結合しているかという、RNA-Seqとは別の情報から選ばれています。
HuR target 200遺伝子は、RNA-Seqではあまり目立たない
この200遺伝子をRNA-Seq上で見てみると、controlとknockoutの間で大きく変化している遺伝子ばかりではありません。散布図では、多くの遺伝子が対角線の近くにあります。HuR targetだけを強調しても、特別な集団としてはほとんど見えません。
一方、200遺伝子をまとめて見ると、knockoutでやや発現が低下する方向への偏りが見えます。著者らは、RIP-Seqで選んだHuR target 200遺伝子と、それ以外の遺伝子について、HuR-FKOとcontrolのfold changeの累積分布（CDF）を作成し、Kolmogorov–Smirnov testで両者を比較しています。その結果、HuR targetは他の遺伝子よりもdownregulated側に偏っていました。個々の遺伝子で大きな発現差を探すのではなく、200遺伝子全体の分布が少しずれているかを見ています。Kolmogorov–Smirnov testは、2群の累積分布の形や位置の違いを評価する検定です。この場合、各遺伝子を個別に有意・非有意に分けるのではなく、HuR targetという遺伝子群全体がどちらの方向へ動いているかを評価しています。
HuR target 200遺伝子のRNA abundance。 左： eWATにおけるWTとHuR knockoutのRNA-Seqデータ。黒点はRIP-Seqから選ばれたHuR target 200遺伝子で、多くは対角線付近にあり、散布図では大きな違いは見えません。右： 各サンプルの発現量分布。
INSIG1はRNA-Seqでは目立つ遺伝子ではなかった
HuR targetの中から、著者らはINSIG1に注目しています。INSIG1は脂肪細胞分化との関係がすでに知られていたため、機能的にも興味深い候補でした。
しかし、RNA-Seqだけを見るとINSIG1の変化はそれほど大きくありません。今回このデータを再解析すると、BAT、eWAT、iWATの3組織で変化の方向はそろっていましたが、t検定ではP<0.05にはなりませんでした。
もしRNA-SeqのDEG解析だけを使い、P値の小さい遺伝子から候補を選んでいたら、INSIG1は候補に入らなかった可能性があります。
eWATにおけるWTとHuR knockoutのVolcano plot。 黒点はRIP-Seqから選ばれたHuR target 200遺伝子、赤茶色はINSIG1を示します。HuR targetの多くはRNA-SeqのDEG解析では特に目立つ位置にはなく、INSIG1も大きな発現差や小さなP値を示す遺伝子ではありませんでした。
INSIG1は、その後の実験で検証された
著者らはINSIG1について、さらに実験を行っています。HuRがINSIG1 mRNAに結合すること、INSIG1の3′UTRに作用すること、HuR欠損によってINSIG1 mRNAの安定性が低下すること、さらにタンパク質量や脂肪細胞分化にも影響することを調べています。
その結果から、HuRがINSIG1 mRNAを安定化し、それを介してadipogenesisを抑えるという機構が示されています。INSIG1は、RNA-Seqで最も大きく変化した遺伝子でも、最も小さいP値を持つ遺伝子でもありませんでしたが、この研究では重要な遺伝子として実験的に検証されています。
Translation Efficiencyにも大きな変化は見えない
この研究ではRibo-Seqも行われており、translation efficiencyについても調べられています。著者らは、HuR knockoutによってtranslation efficiencyに有意な変化は認められなかったとしています。
散布図で見ると、controlとknockoutはよく似ており、HuR target 200遺伝子だけを強調しても明瞭な違いは見えません。今回、同じ200遺伝子を分布で見てみると、わずかな偏りがあるようにも見えます。
ここはRNA-Seqの結果と似ています。RNA abundanceについても、散布図で見れば差はごく小さく、200遺伝子全体の分布を比較することで偏りが検出されました。一方、translation efficiencyについて著者らは、有意な変化は認められなかったとしています。
HuR target 200遺伝子のtranslation efficiency。 左： WTとHuR knockoutのtranslation efficiencyの比較。黒点はRIP-Seqから選ばれたHuR target 200遺伝子で、散布図では明瞭な分離は見えません。右： WTとKOの分布。著者らはtranslation efficiencyに有意な変化は認められなかったとしていますが、分布にはわずかな違いも見えます。
RNA abundanceでは遺伝子群全体の小さな偏りが統計的に検出され、translation efficiencyでは著者らは有意な変化とは判定しませんでした。ただ、P<0.05を境に、前者では現象があり、後者では何も起きていなかったと結論するのは早すぎるかもしれません。
小さな変化が、別の段階で重なることもある
タンパク質量はmRNA量だけで決まるわけではありません。translation efficiencyやprotein stabilityなど、複数の段階の影響を受けます。各段階の変化が小さくても、同じ方向に重なれば、最終的にはより大きな差になることがあります。反対方向に動けば、互いに補償し合うこともあります。
さらに、そのタンパク質が上流の制御因子であれば、小さな変化が多数の下流分子へ伝わり、最終的にはより大きな表現型として現れることもあります。
システム生物学的に考えれば、原因側の変化が小さく、結果側の変化の方が大きく見えることはそれほど不思議ではありません。転写因子、RNA-binding protein、kinase、receptorなどでは、少しの変化が多数の下流分子に伝わります。
RNA-Seqでfold changeが小さいから重要ではない、translation efficiencyの差が小さいから影響はない、と個別の結果だけから判断するのは難しいところです。
200個の中には、ほかにも重要な候補があるかもしれない
この研究ではINSIG1が詳しく検証されていますが、RIP-Seqで得られた200個のHuR targetの中で、INSIG1だけが特別に大きなRNA-Seqの変化を示していたわけではありません。既知の脂肪細胞分化との関係も含めて、INSIG1を詳しく調べる理由があったと考えられます。
200遺伝子全体でRNA abundanceに弱い偏りがあり、translation efficiencyにも小さな偏りが見えるのであれば、この中にはINSIG1と同程度、あるいは別の制御段階ではそれ以上の影響を持つ遺伝子が残っている可能性もあります。
たとえば、RNA量は少し低下するがtranslation efficiencyはほとんど変わらない遺伝子もあれば、RNA量とtranslation efficiencyの両方が少しずつ同じ方向に変わる遺伝子もあるでしょう。RNA量の変化がtranslation efficiencyによって補償される遺伝子も考えられます。
どの候補を次に実験するかは、P値の順位だけでは決まりません。既知の機能や直接結合、RNA abundance、translation、protein levelなどを組み合わせて考える必要があります。
Enrichment解析で見えているのは何か
GSEAは、個々の遺伝子の変化が小さくても、同じpathwayに属する遺伝子がまとまって同じ方向に動けば、その傾向を検出できます。
今回の研究では、RIP-Seqで得られたHuR target 200遺伝子についてGO enrichment解析を行うと、RNA metabolismやRNA processingに関連するprocessが有意に濃縮されました。一方、RNA-SeqのGSEAでHuR knockoutにより上昇したpathwayに含まれる遺伝子と、HuR target 200遺伝子との重なりはほとんどありませんでした。著者らはこの結果から、RNA-SeqのGSEAで見えた炎症関連pathwayなどは、HuRの直接標的そのものではなく、HuR欠損後に生じた下流の変化である可能性が高いと考えています。
Enrichment解析は、どの遺伝子群に対して行うかによって解釈が変わります。一般的なRNA-Seqで、発現差があるとして抽出された遺伝子群に適用すれば、実験条件によって生じた下流の機能変化を捉えられる可能性があります。一方、RIP-SeqやChIP-Seqなどから得られた、原因に近いと考えられる遺伝子候補に適用すれば、その制御因子がどのような分子機構や機能に関わっているかを考える手掛かりになります。
以前、ChIP-Atlasを使って転写因子の結合候補とRNA-Seqの発現変化を比較した解析では、結合候補遺伝子の発現変化は上昇・低下が混在し、散布図や個々のfold changeを見るだけでは明確な傾向を捉えにくいという結果になりました。しかし、今回のHuRの例では、個々の遺伝子の変化は小さくても、HuRの結合候補群とそれ以外の遺伝子群について発現変化の累積分布を比較することで、弱い偏りが検出されています。ChIP-Atlasを使った解析についても、結合候補遺伝子を一つずつ評価するだけでなく、非結合遺伝子との累積分布の違いをCDFで可視化し、Kolmogorov–Smirnov testで比較すれば、これまで捉えられなかった傾向が見えてくるかもしれません。
大事なのは、何を捉え、どうつなぐか
この研究では、RNA-Seqだけで結論を出していません。RIP-SeqではHuRがどのRNAに結合しているか、RNA-SeqではHuRを欠損させた後に発現量がどう変わるか、Ribo-Seqではtranslation efficiencyがどう変わるかを見ています。さらに、protein levelやfunctional assayを使って、その変化が分子機能や表現型につながっているかを調べています。
それぞれの実験で測っているものは違います。研究では、どの現象をどの実験で捉えるかを考え、その結果をどうつなげるかが重要になります。これはデータ解析でも同じです。何を見たいのかが決まっていなければ、DEG解析、GSEA、クラスタリングを順番に実行しても、単なる解析手順になってしまいます。
近年のAIの発展により、データ解析の作業自体は以前より簡単になりました。RやPythonのコードを書くことも、DEG解析やGSEAを実行することも、候補遺伝子を整理したり関連論文を探したりすることも、以前より容易になっています。
一方で、どの実験で何を測るのか、それぞれの結果をどうつなげるのか、どの結果が原因に近く、どの結果が下流なのか、次に何を確かめるべきなのかは、解析を自動化しただけでは決まりません。何を測るかを決め、異なる実験結果をつなぎ、次の問いを作る。研究そのものをどう組み立てるかは、今後も研究者が担う役割であり続けるでしょう。
