【R】遺伝子発現解析を学んでみたい-4③(DESeq2)
・お題:遺伝子発現の多寡を解析するパッケージとして、EdgeR、DESeq2、limmaなどいろいろあるらしい。今回もDESeq2をやってみたい。
・私は素人なので、tutorialをなぞる感じで学んでいきたい。今回のtutorialは以下。正しいことは元サイトを見て頂きたい。
・前回いろいろ描いてみた。今回はエンリッチメント解析をしてみたい。
・チュートリアルでは、gprofilerといパッケージを使っていた。g:Profilerは、遺伝子リストの特徴付けと操作を行うための公開 Web サーバーらしく、視覚化などもできるシンプルなwebインターフェイスらしい。正しいことは元サイトをご確認いただきたい。
・g:Profilerはwebインターフェイスであり、そのR版がgprofilerらしい。調べてみると、バージョンアップしてパッケージ名がgprofiler2になっていた。インストールは割愛。
使い方は以下より。gprofilerと比べると、関数がや機能が結構変わっている。
・まずは、発現誘導がかかっていた遺伝子に着目する。resultDESeq2 から、検索用の遺伝子名のセットを作成する。
>library(gprofiler2)
> resultDESeq2.df <- na.omit( data.frame( resultDESeq2 ) ) #resultDESeq2をdataframe化して、欠損値を含む遺伝子を除去。
> resultDESeq2.df %>% head()
> induced.sign <- rownames(resultDESeq2.df)[resultDESeq2.df$log2FoldChange >= 2 & resultDESeq2.df$padj < 0.01]#log2FC >=2とp-adj< 0.01の遺伝子名を抽出
> induced.sign %>% head()
[1] "yal061w" "yar009c" "yar071w" "ybl005w-b" "ybr012w-b" "ybr067c"
・gostを使ってエンリッチメント解析をやってみる。
> term.induced <- gost(query=induced.sign, organism="scerevisiae")
> term.induced$result %>% head()
・ちょっと図示してみる。
> result_sorted <- term.induced$result[order(term.induced$result$p_value),]
> p <- gostplot(term.induced , capped = FALSE, interactive = FALSE)
> pp<- publish_gostplot(p,
+ highlight_terms = result_sorted$term_id[1:10])
> pp
・term_nameが付くとイメージしやすくなる気がする。
・同様に、発現が下がった遺伝子のエンリッチメントも調べてみる。といっても、遺伝子リスト作成時のlog2FCを>=2から<=-2に変えるだけ。
> repressed.sign <- rownames(resultDESeq2.df)[resultDESeq2.df$log2FoldChange <= -2 & resultDESeq2.df$padj < 0.01]
>
> term.repressed <- gost(query=repressed.sign,
+ organism="scerevisiae",
+ ordered_query = TRUE)
> term.repressed$result %>% head()
> result_sorted2 <- term.repressed$result[order(term.repressed$result$p_value),]
>
> p2 <- gostplot(term.repressed, capped = FALSE, interactive = FALSE)
> pp2<- publish_gostplot(p2,
+ highlight_terms = result_sorted2$term_id[1:10])
> pp2
・工夫が必要かもしれない……。
おわり。