edgeRによる遺伝子発現解析・GSEAとclusterProfierによるパスウェイ解析

日付；2025/02/11（火）

はじめに

先日、RNA-seqの解析でのFASTQのQC、アダプター除去（アダプタートリミング）、ゲノムへのマッピングの方法を記した。ここでは、それに引き続く解析であるedgeRによる遺伝子発現解析の方法、GSEAを使ったGSEA（Gene Set Enrichment Analysis）、clusterProfilerを使ったOverrepresentation Analysisの方法を記すことにする。自分は、ここで記す一連の解析がRNA-seqのデータの解析の基本であると考えている。もしこのRNA-seq解析の目的が、唯のデモンストレーション、例えば、あるタンパク質のノックダウンの効果やそのパスウェイなどが既に他の実験である程度証明出来ており、それを更にRNA-seqで証明するだけ、のような場合、ここに記す基本的な解析で十分な気がする。時間があればそれ以上深堀しても良いのだが、ただ複雑になって説明しにくくなるだけな気がする。

上記とは逆に、探索的なことを行うことで候補となるメカニズムの当たりをつけに行くような場合は上記のような遺伝子発現解析とパスウェイ解析だけでは足りない可能性がある。この場合、自分であればPCA（Principal Component Analysis）やクラスタリングによる候補遺伝子の抽出、転写因子ネットワークの解析、Immune deconvolutionによる腫瘍内免疫環境の推定などを行う。

なお、GSEAに関しては過去に書いているので、それらだけが知りたい場合はそちらの方が有用かもしれない。そこではBroad InstituteのGSEAを使った解析について書いている。他にもシングルサンプルのGSEA（pre-ranked）やRを使ったsingle sample GSEAを書いている。

今回の解析の目的

今回の解析はヒトではなくマウスである。マウスのアノテーションや情報はヒトよりも随分と少ないので、この点はちょっと気をつける必要がある。

今回の解析はシンジェニックモデルから得たいくつかのマウス腫瘍とその宿主（マウス）の全血のRNA-seqであり、その担がんマウスにとある薬剤（CPD）を投与し、そのコントロールであるDMSOと比較する、という実験である。このあたりの薬剤は都合上詳しく言えないが、ここにある方法を参考にする場合は、該当する部分を自分たちの関心のある薬剤に置き換えれば良い。

そして、今回の目的は、とあるパスウェイが薬剤投与群でちゃんと動いているかどうか説明することである。もはや、説明ができ次第、それで解析終了と言える。上記でも述べた通り、このような目的であれば、深追いの探究的な研究は必要ないと言える。深追いの解析をしていると色々な発見があって面白いのだが、そこは時間とか、他のやることとか、コンピューターリソースとか、ライセンスとか、更には費用とはも考えて、色々なバランスで行うか行わないかを決めなければならない。

edgeRとDESeq2のどちらを使うか

edgeRもDESeq2も似たような計算方法を使っているが、自分は主にedgeRを使っている。最終的に個人の好みでDESeq2とedgeRを使い分ければ良いと思っている。

DESeq2とedgeRの主な差別化点としては、DESeq2はきっちり2群の比較しか出来きないが、edgeRの方はより柔軟に解析することができるところと考えている。edgeRは2群の比較にしてもいくつか選択肢（たしかDESeq2もあったと思う）があり、さらに、Anova-likeな検定で複数群の比較も可能である。これはDESeq2にはない利点である。Anova-likeな検定では、複数群で発現に差がある遺伝子はわかるが、どの群で違っているのかはわからない点には注意が必要である。それに加えて、言うても厳密にアルゴリズムを合わせて計算したことは実はないし、詳しい理由を説明する力が自分にはないので確かな保証は出来ないのだが、自分の経験上edgeRの方がDESeq2よりも多くのDEGを拾うことができるように思う。このあたりはおそらく計算に依って異なるのだろうから、アルゴリズムを合わせて計算してみたら、逆にDESeq2のほうが多かった、とかいう場合もあるだろうと思う。その場合は、そうなんだろう。RNA-seqの解析は、大量のデータから候補を拾ってくるあたり、ある意味スクリーニングのような意味合いもあるので、候補遺伝子は多い方が良く、そこから関連のあるパスウェイだって多く拾ってくることが出来るはずである。依って、自分は経験上多くのDEGを拾うことができるedgeRを使うことにしている。以前はDESeq2を使っていたのだが、最近ではどうもedgeRのほうが目的にRNA-seqの目的に合っているように思っている。上記の理由から、様々な状況に対応できるためである。「あの解析ではDESeq2を使ったけど、この解析ではedgeRを使いました」という場合があったが、その時「それやったら前の解析もedgeRでやったほうが良かったかもな….」と思っていた。だったら、はじめから柔軟な解析が可能なedgeRである。

Gene set enrichment analysisとOverrepresentation analysisの違い

パスウェイ解析として主にGSEA（Gene Set Enrichment Analysis）とOverrepresentation Analysisがある。前者は、解析したい遺伝子群（対象となる遺伝子・クエリとして入力する遺伝子群。例えば、2群を比較して統計的有意に発現に差があった遺伝子群）を発現の高い順から低い順に並べて、解析したい遺伝子の発現量が、事前に用意したあるパスウェイに属する遺伝子群（Gene Set）に対しどのくらい寄与しているか、ということを種々の方法により計算し、どちらの群（たとえば薬剤投与群、またはDMSO群）で、どのパスウェイ（Gene set）がより動いているか、働いているか、機能しているか等を評価する方法である。

一方、Overrepresentation AnalysisはGSEAよりも直感的でシンプルな解析であり、解析したい遺伝子群が、どのパスウェイにどれくらい所属しているかを解析する方法を採る。解析結果を見ればなんとなくわかるが、Overrepresentationという呼称は、解析したい遺伝子群が、例えば、GOやKEGGなどに登録されたパスウェイを構成する遺伝子のうちどれくらいを占め、それが他のパスウェイに比べてどれだけ確かかを検定しているところから来ているのだろう。主観的なイメージではあるが、Overrepresentationとは、その言葉通り、「多く（over）に代表とされる（represetation）」パスウェイを見つけるための解析といった感じである。こういった解析なので、Overrepresentation Analysisでは、各遺伝子の発現量やランクは特に必要ない（あってもいいが、それはもうGSEAみたいな計算になるはず。ここで使用するclusterPrifilerは発現量をノード（点・各遺伝子）に色を付けるために使用する。）。一方、GSEAは発現量の値が必要である。この当たりが、GSEAとOverrepresentation Analysisの違いである。

計算方法の違いは、すなわち、使い方の違いでもある。GSEAは対象となる遺伝子の発現量の値を使ってあるパスウェイへの寄与を調べてるので、Overrepresentation analysisよりも解析するための条件が多い（発現量が必要である）。発現量がすごく低かったり、解析したい遺伝子群が少なすぎたりすると上手いこと計算できないし、どちらかと言えば、Overrepresentation Analysisよりも厳しい解析になる用に思う。特に、ここで使用するBroad InstituteのGSEAというソフトは、さらにサンプル間でPermutation testを行うという徹底ぶりである。逆に言えば、GSEAでヒットしてきたパスウェイと、他の実験データとの一貫性があれば、それはもうそのパスウェイの関連性が非常に高いものと考えることが出来ると思う。一方、Overrepresentation Analysisは、各遺伝子がどのパスウェイに所属しているのか程度しかわからないように思う。だから、GSEAに比べてヒットしてくるパスウェイの数は圧倒的に多い。しかしながら、たとえば、細胞Aでは細胞Bに比べて本当に現象Xに関するパスウェイが深く関わっているかどうかを単に示すだけ、のような解析ではOverrepresentation Analysisでも良くて、むしろ直感的に分かりやすくて良いと考えている。その理由は、Overrepresentation Analysisは、比較群に対して何らかの顕著な違いをもつ遺伝子、例えば、RNA-seqであればDEGを使用するため（これはGSEAでもそうであるが。そうしないと絞り込みが非効率である。）である。逆に、そういったフィルターなしにOverrepsesentation Analysisを使う意味はほとんどない。Research Misconductにならないようならば、この辺は上手く使っていけば良いと考えている。

前置きがかなり長くなってしまったが、以下からRを使って解析を進めていく。Rはかなり会話的なソフトなので、各ステップ毎に何をやったかを説明していく。