研究に役立つ！enrichplotでGO解析を視覚化する方法【其の弐】

バイオインフォマティクス 統計解析ソフトR エンリッチメント解析

アメリエフ解析チームのkutsuyです。
前回に引き続き、enrichplot を用いたGO解析結果の可視化について検討してみます。

enrichplotでGO解析を視覚化する方法【其の壱】
enrichplotでGO解析を視覚化する方法【其の弐】　←本記事

アメリエフのRNAseq受託解析ではエンリッチメント解析までを標準の納品物としております。 興味のある遺伝子リストでエンリッチしている機能から新しくわかる発見もあります。 文章も長くなりますので、リラックスしてご覧ください。 以下のサイトを引用しております。

yulab-smu.top

動作環境

  • R version 4.3.0
    • DOSE v3.26.2
    • enrichplot v1.20.3
    • dplyr v1.1.4

Dot plot の出力

前回はドットプロットを作成しましたが、mutate関数を使用した以下のコマンドも実行してみます。

> mutate(edo, qscore = -log(p.adjust, base=10)) %>% barplot(x="qscore")

ここで手元の環境では下記のようなエラーが出ました。

Error in mutate(edo, qscore = -log(p.adjust, base = 10)) %>% barplot(x = "qscore") : 
  could not find function "%>%"

これはエラーの内容からもわかるように、「%>%」が認識されていないためです。 %>%はパイプ演算子とも言われ、%>%の左側のコマンドの出力結果を右側のコマンドに引き渡す機能になります。 補足ですが、R界隈で有名なtydyverseというパッケージにdplyrが含まれているため、tydyverseを読み込むことでもパイプ演算子が有効となります。

> library(dplyr)

今回はパイプ演算子を有効にするために、dplyrパッケージを読み込みました。 再度、mutate関数のコマンドを実行してみたいと思います。

> mutate(edo, qscore = -log(p.adjust, base=10)) %>% barplot(x="qscore")

ここで手元の環境では下記のようなエラーが出ました。

Error in UseMethod("mutate") : 
  no applicable method for 'mutate' applied to an object of class "enrichResult"

edo というオブジェクトは enrichResult オブジェクトという特殊な構造になっています。 mutate 関数は第一引数にデータフレーム型のオブジェクトを設定する必要があります。 では、as.data.frame 関数でデータフレーム型に変換してみるとどうでしょうか。

> mutate(as.data.frame(edo), qscore = -log(p.adjust, base=10)) %>% barplot(x="qscore")

ここで手元の環境では下記のようなエラーが出ました。

Error in barplot.default(., x = "qscore") : 
  argument 2 matches multiple formal arguments

mutate関数は実行されているのですが、その出力を barplot 関数の入力として実行する際に、 x = "qscore" という引数が複数の引数にマッチしてしまっているためエラーとなっています。 また、ここで実行されている barplot 関数は base パッケージ（デフォルトで使用可）であり、 enrichplot パッケージの barplot 関数ではないことにも注意が必要です。

データを整形して ggplot2 などのパッケージで同様の作図を試みることは可能ですが、かなり話が長くなっているので割愛します。
ここで時間をかけてまで私がお伝えしたかったことは、「ネット上のコマンドは再現できないことがある」ということです。
パッケージのバージョンが揃っていないことが原因なこともありますし、コマンド自体が間違っている場合もあります。 今はAIが主流になっているので、エラーを解決するよりも代替案をAIに聞いた方が早い場合が多々あります。 但しエラー解決に時間をかけて取り組むと知識が上達するメリットもありますのでいろいろチャレンジしてみてください。

さて、作図が全然進まないので次の項目もどんどん実行していきましょう。
前回の実行の後に一度 R を終了している場合は、library 関数で DOSE パッケージを再度読み込んでから実行してください。library 関数は R 立ち上げのたびに実行する必要があります。

> edo2 <- gseDO(geneList)
> dotplot(edo, showCategory=30) + ggtitle("dotplot for ORA")

RStudio などの総合環境で実行されている方はプロットが表示されたのではないでしょうか。

Dotplot

前回の barplot では p.adjust が小さい順に並んでいましたが、今回は GeneRatio の順序で並んでいますね。 また、前回もでしたが showCategory という引数でタームの数を指定することができます。 enrichplot は ggplot2 ベースのグラフで可視化するツールですので、ggplot2 の構文でカスタマイズできます（色、形など）。 ぜひ自分のデータを分かりやすいグラフに加工してみてください。

Networkplotの出力

次は、Networkplot を作図してみましょう。 まず、setReadable 関数で edo オブジェクトの遺伝子ID（ENTREZID）を遺伝子シンボルに変換します。

> edox <- setReadable(edo, 'org.Hs.eg.db', 'ENTREZID')

エラーなく実行出来たら、次に cnetplot 関数を実行して、その結果を p1 という新しいオブジェクトに代入します。

> p1 <- cnetplot(edox, foldChange=geneList)
> p1

ここで手元の環境では下記のような警告（エラーではなく）が出ました。

Warning message:
In cnetplot.enrichResult(x, ...) :
  Use 'color.params = list(foldChange = your_value)' instead of 'foldChange'.
 The foldChange parameter will be removed in the next version.

これは Warning（警告）なのでコマンドの実行自体には影響はありません。 foldChange という引数は今後パッケージが更新されると無くなるので注意喚起をしてくれています。 p1 というオブジェクトには下記のプロットが含まれていました。

Network plot

ここでお気づきの方がいらっしゃるかもしれませんが、WEB上に掲載されている図とやや違う気がします。 作図の一部、特に回転系の図については同じ意味を持つ図でも違う形になることがあるので注意してください。 研究者の方は図表の再現性に強い関心があるかと思います。シード固定などで再現性を担保できることがあるのでぜひ検索してみてください。

サイト上にもありますが、シード固定の例は以下のようなコマンドになります。 カッコの中の数字は任意の数字を入れてください（123である必要はありません）。 作図が多いので割愛しますので、ぜひいろいろ実行してみていただければと思います。

> set.seed(123)

Heatmap plot の作図

次は、Heatmap plot を作図してみましょう。 こちらはそのままのコマンドで行けるはずです。

> p1 <- heatplot(edox, showCategory=5)
> p2 <- heatplot(edox, foldChange=geneList, showCategory=5)
> cowplot::plot_grid(p1, p2, ncol=1, labels=LETTERS[1:2])

Heatmap plot

遺伝子とタームの関係性がとてもわかりやすいですね。 着目すべき遺伝子のピックアップなどに利用できそうです。

Tree plot の作図

次は、Tree plot を作図してみましょう。 こちらもそのままのコマンドで実行できることを確認しました。

edox2 <- pairwise_termsim(edox)
p1 <- treeplot(edox2)
p2 <- treeplot(edox2, hclust_method = "average")
aplot::plot_list(p1, p2, tag_levels='A')

Tree plot

ターム（機能）がクラスタリングされています。 二群間で有意に変化している機能が多いクラスタが一目でわかるような図ですね。

内容も長くなってきましたので、次の作図で今回は一区切りにしたいと思います。 最後までどうぞお付き合いください。

Enrichmentmap plot の作図

次は、Enrichmentmap plot を作図してみましょう。 図が多いので一つだけ作図してみます。 言い遅れましたが、これまでは複数の図を別々のオブジェクトに代入し、 cowplot パッケージの plot_grid 関数で並べて表示していました。

> edo <- pairwise_termsim(edo)
> p1 <- emapplot(edo)
> p1

Enrichment map plot

二群間で有意に変動していたタームが一つ一つのプロットで表され、 含まれる遺伝子が共通すると線で結ばれる仕様となっています。 よって、似たような機能は蜘蛛の巣のようにクラスタとなるため、 大まかにどのような機能分類があるかわかりやすいですね。
くどいようですが、シード固定をしなければ実行のたびに図は回転します。
clusterProfiler というこの分野で有名なパッケージがあるのですが、 そちらを使用したコマンドも紹介されていますので是非トライしてみてください。

まとめ

今回は前回に引き続き、enrichplot を用いた作図を試してきました。
それぞれの図には意味がありますので、作図が出来たことに満足せずに（私はすぐに満足しがちですが）、 図から読み取れる生物学的な意味を確認し、考察を加えることが非常に重要です。
というか考察するために作図をしていますので、作図が目的とならないようにご留意ください。 作図はあくまでも手段です、目的はそれらのデータから新しい知見を得ることだと思います。

ということを目的を忘れがちな今の私に一番言いたいです。
次回でこのシリーズは終わる見込みですので、引き続きお付き合いくださいませ。