いろいろ倉庫

KNIME、EXCEL、R、Pythonなどの備忘録

【R】遺伝子発現解析を学んでみたい-4③(DESeq2)

・お題:遺伝子発現の多寡を解析するパッケージとして、EdgeR、DESeq2、limmaなどいろいろあるらしい。今回もDESeq2をやってみたい。


・私は素人なので、tutorialをなぞる感じで学んでいきたい。今回のtutorialは以下。正しいことは元サイトを見て頂きたい。

dputhier.github.io

・前回いろいろ描いてみた。今回はエンリッチメント解析をしてみたい。

チュートリアルでは、gprofilerといパッケージを使っていた。g:Profilerは、遺伝子リストの特徴付けと操作を行うための公開 Web サーバーらしく、視覚化などもできるシンプルなwebインターフェイスらしい。正しいことは元サイトをご確認いただきたい。

biit.cs.ut.ee

・g:Profilerはwebインターフェイスであり、そのR版がgprofilerらしい。調べてみると、バージョンアップしてパッケージ名がgprofiler2になっていた。インストールは割愛。

biit.cs.ut.ee

使い方は以下より。gprofilerと比べると、関数がや機能が結構変わっている。

cran.r-project.org

・まずは、発現誘導がかかっていた遺伝子に着目する。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

・工夫が必要かもしれない……。

 

おわり。