RNAシークエンス解析関連のおすすめライブラリ

日付；2024/05/19（日）

何かの標的分子と疾患との関連を知りたい場合、結局のところ通常のRNA-seqなどを利用する場合が多い。なので、一般的な遺伝子発現解析ツールも準備する必要がある。最近ではRNA-seqもかなり普通になってきたので、以下のライブラリで通常のことは出来るように思う。ちなみに、以下のソフトを使った解析で一貫性のある結果が得られた場合、それは生物学的に有意であり、すぐにでも生物実験にトランスレーションできるような結果であるのではないかと思う。逆に、非常に曖昧で、統計的に緩い条件にしなければならないような結果だった場合、解析による深追いは無駄になるように思うし、もしかしたら解析の条件が間違っているか良くないとか、条件設定が悪いとか、技術的な問題があるとか、そういった可能性がある。

正直、通常のRNA-seqの高次解析についてはRが良いと思う。もはや解析手順が統計的に確立されているし、高次解析で時間がかかるのは後述するbaySeqくらいである。それ以外のFastqファイルのQCとかプロセッシングは、Rではなくシェルである。Rへの実装もあるようだが、シェルを使って行われているのが一般的である。ということで、個人的にRNAシークエンスによる遺伝子発現解析に使用しているライブラリである。以前、ここにけっこうちゃんと書いた。

QC・Preprocessing

多くの場合、FASTQファイルが返ってくると思う。その場合、まず最初にリードのQC（Quality Check）を行い、必要に応じたトリミングを行うのが一般的である。これは、リードのPreprocessingとか呼ばれていると思う。おそらく、世の中にはたくさんのパッケージやライブラリがあるだろうが、以下の表にあるソフトウェアで大体のことは出来るはずである。

シークエンスをしてくれる業者によっては、もはやQCやProcessingが必要ないくらい綺麗なFASTQファイルが返ってくることがある。FastQCをやってみて明らかに綺麗なリードである場合、無理にトリミングなどはする必要がないと思う。一方、SRA（Sequence Read Archives）などからダウンロードしたデータなどでは、ライブラリ調整キットのバージョンが古かったり、そもそもライブラリのクオリティーが低かったりで、この辺りのQCやPreprocessingが必要な事が多いように思う。

あと、後述するマッピング後のBAMファイルのQCを行う場合は、rseqcが有用である。マッピング後のクオリティーがわかる。QCの過程でPCAをやってみたらあるファイルだけ主成分がずれていたとか、あるファイルだけカウントデータがおかしいとか、そういうことが詮索できる。

マッピングとカウント

QC一通り終わったらリファレンスゲノムを使ってリードをマッピングしていくのだが、そのためにはやはり定番のSTARが良いと思う。STARはメモリは爆食いするが最も認められているマッピングツールであると思う。特に、fusionChatcherなどではSTAR2パスなんかが使われているようなので、互換性も良いだろうと思う。もちろん、HISAT2やKallistoなどでも良いと思う。HISAT2とかも、リファレンスゲノムをちゃんと自分で構成してやれば信頼できる解析になるのではないかと思う。

つぎにマッピングされたリードをカウントするが、このためにはfeatureCount一択じゃあないだろうか。他にはsalmonとかKallistoなどが考えられるが、これらはちょっとアルゴリズムが違うように思う。あとは、以前Nature Protocolに出ていように、HISAT-StringTie-Ballgownを使ってマッピングからカウントまで行っても、わたしは一向にかまわん。

RNAシークエンスの解析をやっていると、DEG（Differential Gene Expression）の解析はもちろん、融合遺伝子の解析なんかもやりたくなる。その場合はfusionCatcherが良いだろうと思う。他のパッケージもあるが、ベンチマークが結構良いのがfusionCatcherである。ただし、このパッケージもだんだん古くなってきているので要注意である。これに必要なPythonなんか、 ver. 2.7である。正常に動けば問題はないが、これはつまり開発がもはや止まっている事を示している。

STARとかfusionCatcherをインストールする場合、samtools、BWA、BOWTIE2とかも必要になるはずである。

DEG（Differentially Expressed Gene）解析

リファレンスゲノムへのリードのマッピングが終われば、次はDifferentially Expressed Gene（DEG）を見つけにいく解析である。これも色々なライブラリがあるが、おそらくedgeRかDESeq2が定番だと思う。両者は似たようなNormalizationを行うので、好みに応じて使って良いように思う。詳しいアルゴリズムによるのだろうが、似たような計算なのになぜかDESeq2の方が検定が厳しいような印象がある。なので、スクリーニングと割り切って多くの候補遺伝子が必要な場合なんかは、おそらくedgeRの方が良いように思う。どうせqPCRで拾ってきた遺伝子の発現量をバリデーションしなければならないので、そこは多くの候補を拾ってきた方がお得である。もう一つの大きな違いは、edgeRではAnova-likeな解析を行う事ができる点である。例えば、グループが3群以上あった場合に、各々の遺伝子発現量がどこかの群で違う可能性を示してくれる検定方法である。このように考えると、RNAシークエンスのDEG解析にはEdgeRを使う方が汎用性があるように思う。だから、自分はよくedgeRを使っている。ちなみに、3群以上ある場合は、得られたp値をp.adjust()で補正する必要がある。これをやってしまうと、結構厳しい検定になってしまうし、RNAシークエンスはスクリーニングという意味合いが強く、後にバリデーションする必要があるので、3群以上あった場合でも関心のある群とコントロール群の2群の比較をやってしまっても良いように思う。正直にp値を補正して、全く有意な遺伝子が取れてこないとか、なんか微妙である。

比較群が多く、検定の多重性が問題になってしまう場合は、baySeqを使う。これは多群の比較をしてくれる。最悪、Pairewiseに比較してくれるので、とても良い。ただし、計算時間がかなり長い。最近のハイスループットシークエンサーだとリードの数もかなり多い場合があるが、それだと平気で数時間かかったりする。余談だが、やはり無駄にディープシークエンスは良くない。目的に対して適切なデプスでシークエンスするべきである。計算時間が長いので、そこはどの群が必要なのかちゃんと考えて使用する必要がある。

GTFファイルを読み込んで、後に行うエンリッチメント解析の準備もしたりする必要があるので、それにはrtracklayerを使う。個人的には、GTFファイルはGENCODEのものを使うのが良いと思う。Ensemblとかの方が、汎用性は高いような気もするけど。

エンリッチメント解析

DEG解析が終われば、パスウェイ解析を行う。これはエンリッチメント解析とか、オーバーリプレゼンテーション解析とか言われている。エンリッチメント解析の方が一般的な呼び方であると思うが、けっこうフワッとしているように思う。

そのエンリッチメント解析の中で、個人的には最初にGSEA（Gene Set Enrichment Analysis）を行うことにしている。これは、DEGとして得られた遺伝子のうち、どのくらいの数の遺伝子がどのパスウェイに属しているのか解析する方法である。代表的なRのライブラリとして、まずはclusterProfilerを挙げる事が出来る。これは非常に高機能なライブラリで、GSEAだけでなく、 DAVIDのようなエンリッチメント解析も行う事が出来る。個人的に、これが普及している理由の一つが、GOだけでなく、KEGGやReactomeなども利用することが出来る点だろうと思う。特に、KEGGとpathviewも利用できるところがかなり強い。しかし、なぜか動作が不安定で、サーバーに繋がらなかったりする事が多いので、正直、シームレスではないように思う。だから最近では第一選択では無い。他にパッと思いつくライブラリとしてはfgseaである。どちらかといえば、fgseaが使い方もシンプルで好きである。

エンリッチメント解析で圧倒的に信頼できるのが、Broad InstituteのGSEAである。上記のソフトと比較して、メンテナンスもされているし、統計もしっかりしている印象がある。一番の違いは、このGSEAは、 ノーマライズされたカウントデータを直接利用し、必要に応じてpermutation testも行う事が出来る点である。clusterProfilerもfgseaも、カウントデータを用いたエンリッチメント解析ではなく、DEGのLog2FC値や、ランキングの値だったりする。GSEAにもこのランクを使った解析方法（Preranked GSEA）も出来るが、これをやってみると感度が高すぎて、物凄い数のパスウェイがエンリッチされてくるので、要注意である。あんなにたくさん統計的有意とされるパスウェイが拾われてくると、解析がかなり大変だし、こっちとしてはある程度候補遺伝子を絞り込むためにこういった解析をしているのに、その意味が薄れるような気がする。あと、このGSEAの図は、論文でけっこう見かける。そういう意味でも、Broad InstituteのGSEAが第一選択である。Broad InstituteのGSEAの結果の可視化にはCytoscapeと、そのプラグインであるEntichmentMapが必要である。これは、GSEA自体がうまく行っていればボタンをクリックするだけで生成される。この結果から、任意のパスウェイにおけるLeading edgeを拾ってきたり出来る。Leading edgeというのは、パスウェイのエンリッチメントへの寄与が大きい遺伝子のことである。あるDEGがたくさんのパスウェイに属しているならば、その遺伝子がleading edge geneということになる。こうやってさらにエンリッチしてくことで、ある要因が影響を与えるパスウェイを特定していくことができる。

エンリッチメント解析を行うにあたっては、遺伝子のAccession numberをGene symbolを変換したり、その逆にしたりとする場合がある。なので、必要に応じてbiomaRtやclusterProfilerのbitr()が必要になるかもしれない。もしかしたらnitchnetrでヒトの遺伝子名をマウスのオルソログに直したり、その逆をやったりといった作業も必要になるかもしれない。

個人的に非常に強力だと思う解析がARACNEとVIPERによる影響の大きい転写因子やシグナリングのHubとして機能している遺伝子の推定である。もしDEGで得られた遺伝子や遺伝子群が転写因子が関連するネットワークであれば、インパクトの高い解析になると思う。ただし、残念ながらアカデミックでの利用に限られる。

可視化

DEGやleading edge geneの可視化にはヒートマップが最強である。特に、個人的にはComplexHeatmapというライブラリが欲しい図を作成できる。pheatmapの方が使い方がシンプルで必要最低限ができるので、単純なヒートマップはpheatmapが一番良いと思う。しかし、ヒートマップの上や左右に複数の群のラベルを表示するような場合は、このComplexHeatmap一択である。

また、エンリッチメント解析の可視化には、結局ggplot2を使う羽目になる。clusterProfilerで解析から可視化まで一貫して行っても良いのだが、clusterProfilerがシームレスに欠けるので、図に使うデータを取得できる以上、再現性が良い図を書くにはggplot2が良いと思う。

まとめ

このあたりの解析は今や定番になってきて、すごく普遍的なものになってきてる感じがする。一方で、TCGAなどにはかなりの数の検体数が登録されていることもあり、基礎データのトランスレーションなどではとても有用な解析になる。なので、しっかりと理解しておいてすぐに使えるようにしておいたほうが良いように思う。