日付；2022/10/22（土）

先日、Single Sample（N=1）のRNA-seqは意味があるのか無駄なのかを書いたが、そうは言っても、いろいろな事情によりそのような解析が必要な場合もある。それに、うまく行けばのサンプルの遺伝子発現プロファイルからそのサンプルがどのような細胞や組織に属する推測できるかもしれない。「最終的にお蔵入りになるから無駄だ」とか言っている自分も、結局のところその解析が必要だった。それに、この手の解析、特にGSEA（Gene Set Enrichment Analysis）の手順なんか、いちいち忘れてしまう愚かしい頭脳の持ち主である。正直こういうことを浅はかな知識で公開するのは良いことなのかわからないが、覚書くらいにはなるだろうと思ったので記すことにする。これはSinlge Sample（N=1）でGSEAをやる方法である。ちなみに、既にソフトのインストールは済んでいるものとして記述する。インストールは簡単なので大丈夫だろう。さすがに。

GSEAのマニュアル

マニュアルはここ。各種ファイルフォーマットについてはここ。変な解析をやってしまう前に、一度これは読んだ方が良いと思う。しかし、いちばん重要な結果の見方がしっかり書かれていないように思うのは自分だけだろうか。いずれにせよ、読んでいて損はないと思う。

GSEAのダウンロード

ここからダウンロードする。ダウンロードするためには自分のメールアドレスだのを登録する必要がある。WindowsとMac版などがあるが、個人的に自分のコンピュータに合わせれば良いと思う。コマンドラインもある。Rのバージョンもあるが、２０１９版（２０２２年１０月の時点）らしいので、もしかしたら少し注意が必要な可能性もある。自分はR版は使ったことがないので、よくわからない。というか、２０１９年版のヤツなんて、あんまり使いたくない。

必要なファイル

rnkファイル

Single SamplのRNA-seqの場合はPreranked GSEAを使う。なぜかと言うと、Nが１しかなく、GSEAのソフトではコントロール群と比較したランクが作れないためだ。そのために自分の知恵と知識に従ってランキング付けした遺伝子セットが必要になる。

まず用意するファイルは拡張子が.rnkのrnkファイルである。これは上述の自分の知恵と知識に従ってランキング付けしたRNA-seqのファイルになる。ランキング付けには、色々な方法がある。おそらくTPM（Transpcript per million）で標準化した発現量、コントロール群と比較した発現比（一般的にはLog 2 Fold Change；Log2FCだろう）、さらには、Signed p valueという方法もある。個人的には、N=1の場合は、Log2FC一択だろうと思っている。もちろん、このLog2FCの計算にはTPMにより標準化した遺伝子セットを使う。フォーマットは以下の画像のような感じになる。一行一列目が、#1、そのとなりが空白である。これは最悪エクセルかテキストで作ってしまっても良いと思っている。自分はエクセルを開くのがこの上なく嫌なので、Rで作成してしまっている。

Chipファイル

これは、使用するrnkファイルに書いてある一列目の遺伝子IDなりプローブIDなりプローブ名を、GSEAが理解できる遺伝子名に直すためのファイルである。これは解析のたびにGSEAのソフトがサーバーからダウンロードしてくるので、マウスやヒトなどEntrezやEnsemblの遺伝子IDなどであれば特に用意する必要がない。しかし、自分がマッピングに使ったGTFファイルから遺伝子IDと遺伝子名を抽出してきて使うのが、一番互換性が良く、無駄もない。具体的にはGSEAがダウンロードしてくるファイルは最新ではなく、一部の遺伝子が解析から外れてしまったということを出来る限り防ぐことが出来ると考えている。なので、自分で作る。Rを使って以下のようなファイルを作れば良いと思う。

GMTファイル

これはパスウェイと、それに属する一連の遺伝子セットが書かれたファイルである。これはGSEAやMsigDBからダンロードしてくることが出来る。いちいちダウンロードするのが面倒なので、使用するものは自分のコンピュータに保存したほうが良いと思う。

GSEAする

Load data

上記で作成もしくはダウンロードしたファイルをソフトに読み込ませる。「Load data」から読み込む。問題がなければNo Errorという旨が書かれたウィンドウが出てくる。

Run GSEAPreranked

次に「Run GSEAPreranked」に行く。流れとしては、前のステップでソフトに読み込ませたファイルを選んで、解析の設定をすこしイジって、最後に下の「Run」を押すことである。

Gene sets database

Gene sets databaseには、事前にダウンロードし、前のステップで読み込ませたGMTファイルを選ぶ。ヒトやマウスの場合は、ここでダウンロードできるので、前のステップで読み込ませる必要はない。しかし、それはすごくオープンになっているデータセット、例えばGOとかReactomeとかのデータなので、それ以外にCurationされた遺伝子セットが必要ならば、MsigDBから事前にダウンロードするか、自分で作成して前ステップのように読み込ませる必要がある。自前で作成して事前に読み込ませたGMTファイルは、Gene matrix（local gmx/gmt)のタブのところにある。

Ranked list

次にRanked listから、解析対象となるファイルを選ぶ。これは事前に読み込ませたrnkファイルである。データが問題なく読みこまれているならば、ドロップダウンで出てくる。

Collapse/Remap to Gene symbol

Chipファイルに、rnkファイルに書かれた遺伝子ID（上述の例ではEnsemble gene ID、バージョンなし）をGSEAを行うときに遺伝子名に変換するためのドロップダウンである。デフォルトで良いと思う。ここは気を使う場所であり、rnkファイルとChipファイルの遺伝子IDが対応していない場合、エラーになってしまうので注意である。もちろん、Collapseしないという作戦もある。その場合はrnkファイルにそのまま遺伝子名などを使えば良いのだろうと思う（たぶん）。しかし、RNA-seqの解析ではやっぱり遺伝子IDなんかを使うだろうから（重複しないため）、やっぱりここも必要であろう。これを書いていて思うのだが、やっぱりRNA-seqのマッピングやリードカウントで使ったGTFファイルから引っ張ってくるのが良いように思う。

Chip platform

読み込ませたChipファイルを選ぶ。

Analysis name

次にAnalysis nameで区別できる名前を付ける。これが解析後に作成されるフォルダ名になる。

Save results in this folder

結果は「Save results in this folder」のところで指定したパスに保存される。

Run

設定が終わったら、下にある「Run」を押す。問題なく解析されていれば緑字でSuccessと左下に表示され、何かあればSuccess(with warning)とピンクで示されるが、まぁ、これで良い。Errorとときは「Error！」と赤だった気がする。保存された結果は、「Show results folder」にあるフォルダー（「Save results in this folder」のところで指定したパス）に保存されている。Successのところをクリックすると、結果に飛ぶことが出来る。

GSEAの結果

主な図

よく見る、こんなヤツである。

これは発現比をランキングにしてGSEAを行ったときの結果の一つであり、ランキングの高かった遺伝子、つまり発現が高かった遺伝子群があるパスウェイに属する遺伝子群とどのくらい一致していたかということを示している。残念ながらPreranked GSEAではサンプルの区別がない（N=1だから）ので、”na_pos”(ランキング、この場合は、発現比が高かった遺伝子）と”na_neg”（ランキング、この場合は、発現比が低かった遺伝子）でどのくらいそのパスウェイにマッチしていたか、ということを示す。これは、ランキング（発現比）が高い遺伝子と、データベースにあるパスウェイに属する遺伝子群がマッチしていました、といことを示している。

おそらく結果を読むには注意が必要

しかし、目的のパスウェイがエンリッチされてきたとしても、手放しに喜ぶのは早いように思う。イマイチ詳しくないのだが、結果にはPermutationしたときの解析っぽい（あまり良く理解していないので、こういう状態で公に書くには良くないとは思っているのだが…）図も載っており、メインの図ではエンリッチされているのに、その結果を見ると逆になっていることが良くある。これはおそらく、Permutationをするとエンリッチしない場合もあるってことだと思う（間違っていたら申し訳ない）。やっぱり、メインの図と計算途中の図が一致しているパスウェイが良いと思う。厳しく考えるならば。

Leading edge analysis

これは前述で得られたGSEAの結果のまとめを表してくれるところである。遺伝子を横軸、パスウェイを縦軸にとり、そのクラスタリング結果をヒートマップで示してくれたりするが、正直、GSEAPrerankではヒットする遺伝子セットがかなり多いので、見るのがとても非効率であり、そのためにあまり意味がない。空欄にGSEAの結果が保存されたフォルダーを入力すれば、自動的にまとめてくれる。念の為にGSEAの結果のまとめを保存しておこう、くらいの感じあり、実質このデータを使って何かする、というのはなかなかないと思う。

Enrichment Map Visualization

GSEAの結果得られたパスウェイについて（遺伝子ではない）クラスタリングしてくれ、それをネットワークで示してくれる図である。この機能を利用するにはCytoscapeをインストールし、事前に起動させておく必要がある。しかし、特に起動していなかったとしても、「起動してな」という警告があり、先に進めないだけである。設定としてはP-value Cutoffを0.05、FDR Q value Cutoffを0.05にするのが良いと思う。GSEAはFDRを0.25くらいでアクセプトするとそのマニュアルに書いているが、それだと結構大きいように思う。確かに、FDR Q valueが0.05だと厳しいが、Preranked GSEAはエンリッチメントされてくる遺伝子セットがかなり多いので、厳し目に見るべきだと思う。大体、甘目にみるのはどうなんだろうか。通常のGSEAだって、出来るならば厳し目が良いと思う。

ノードとエッジについての設定は、jaccard+Overlap Combinedのデフォルトで良いと思う。ちなみに、「ノード（Node）」が、丸 、すわなち、その遺伝子セットを示し、「エッジ（Edge）」がそのノードをつなぐ線である。別名は「リンク」でり、リンクのほうが圧倒的に理解しやすい。Edge Computingとか、そういうコンピューティングの分野からのネーミングだろうが、はっきり言ってこれは直感的にもわかりにくいと思う。Edgeはリンクなので、あるノード（遺伝子セットやパスウェイ）に含まれる遺伝子と他のノードに共通した遺伝子がある場合、ノード同士が繋がることになる。

Prerank GSEAで得られたネットワーク（Enrichment Map）をCytoscapeで解析するときの注意

Preranked GSEAの結果得られたネットワークをCytoscapeを使って解析する場合の注意としては、GSEAに使用した遺伝子発現量が反映されていないことである。その解析の結果、何か興味深いネットワークが抽出できたら、そのネットワークに属する遺伝子と実際のデータセットの遺伝子をマッチングして、ちゃんと発現量なり発現比なりを確認する必要がある。つまり、Enrichment Scoreへの寄与が小さい遺伝子もCytoscapeで表示されているパスウェイには表示されている、ということである。Cytoscapeからエクスポートしたテーブルから遺伝子を抽出してきても、もしかしたら発現の変動が小さい遺伝子がほとんどだった、なんてことになりかねない。通常のGSEAの結果得られたネットワークは、そんなことはなかった思う。考えてみればそれは当然で、Preranked GSEAでは、GCTファイル（各サンプルの発現量を記載したファイル）を使っていないので、それはできないってことだろう。

GSEAのrnkファイルを作成するときの注意

注意しなくてはならないのはSigned p valueというヤツを使ってGSEAを行う場合である。個人的に、なぜこんな方法がまかり通っているのか全く理解できない。というか、利用はやめたほうが良いのでは無いかと思う。このSigned p valueというヤツはLog2FCにp値にマイナス１を掛けてかつLog10で変換した値を掛けた値である。これだけ聞いたら、なんか専門っぽい響きではあるが、ちょっと考えてほしい。出来上がった値は、一体どんな値なんだろうか。log2FCだけなら、その単位は発現比であり、それは直接生物実験に使用できしょうなものである。でもそれにp値を正の値をLog10を掛けた時点で、生物学的な意味が消える。この値を使って得られたGSEAに基づいて次の解析をする場合、元のLog2FCを注意深く比較し直さなくてはならない。ということは、Signed p valueを使って得られたGSEA自体に、意味はあったのだろうか。それを使ってなにかのパスウェイが統計的有意にエンリッチされてきたとしても、実際のリードや発現比が低すぎて使えないなんてことになったら、どうするのだろうか。

Preranked GSEAについて思うこと

Preranked GSEAがかなりのパスウェイがヒットしてくるので、解析がかなり大変である。そもそも、発現解析の時点で有意差が出ないってのは、どうしても判断が客観的ではなく、好きになれない。かなりの量のパスウェイがヒットしてくるので、EnrichmentMap（Cytoscape）の結果の解析も、結構たいへんである。