hello，今天小云來(lái)教大家如何對(duì)KEGG富集分析的結(jié)果通過(guò)ggplot2實(shí)現(xiàn)可視化并繪制對(duì)應(yīng)的KEGG通路網(wǎng)絡(luò)圖 ，讓你的分析結(jié)果更加美觀！感興趣的話就和小云一起來(lái)看一下吧！

準(zhǔn)備需要的R包

先來(lái)簡(jiǎn)單介紹一下我們今天用到的R包，分別是DOSE、clusterProfiler、org.Hs.eg.db、ggplot2、tidyverse以及enrichplot包。同學(xué)們要提前下載好哦！在這里小云為大家簡(jiǎn)單的介紹一下其中幾個(gè)R包：

1.?DOSE：DOSE包主要用于基因功能富集分析，它能夠通過(guò)比較輸入的基因列表與已知的基因功能注釋信息，找出在特定功能類(lèi)別中富集的基因。

2.?clusterProfiler：clusterProfiler是一個(gè)用于生物信息學(xué)分析的R包，主要用于功能注釋和富集分析。它提供了針對(duì)多種生物數(shù)據(jù)類(lèi)型的功能富集分析方法，包括GO、KEGG、Reactome等。該包還提供了可視化工具，使用戶能夠更好地理解和解釋功能富集結(jié)果。

3.?org.Hs.eg.db：org.Hs.eg.db是一個(gè)R語(yǔ)言中用于人類(lèi)基因注釋的數(shù)據(jù)庫(kù)包。它提供了對(duì)人類(lèi)基因組的注釋信息，包括基因ID、基因名稱(chēng)、基因符號(hào)、通路信息等。使用org.Hs.eg.db包，用戶可以將基因ID轉(zhuǎn)換為易讀的基因符號(hào)，從而更容易理解和解釋功能富集結(jié)果。

4.?tidyverse：tidyverse是一個(gè)由多個(gè)R包組成的集合，旨在提供一套一致且易于使用的數(shù)據(jù)處理和可視化工具。這些包包括ggplot2、dplyr、tidyr、readr等。其中，ggplot2是用于數(shù)據(jù)可視化的重要組成部分，提供了強(qiáng)大且靈活的繪圖功能，使用戶能夠創(chuàng)建高質(zhì)量的圖表，用于數(shù)據(jù)探索和展示。

5.?enrichplot：enrichplot是一個(gè)用于生物信息學(xué)分析的R包，主要用于功能富集結(jié)果的可視化。它提供了多種繪圖函數(shù)，可以繪制基因集的富集分析結(jié)果，如GO條形圖、KEGG通路網(wǎng)絡(luò)圖等。enrichplot提供了豐富的參數(shù)設(shè)置，以滿足用戶對(duì)結(jié)果可視化的個(gè)性化需求。

KEGG富集分析

enrichKEGG進(jìn)行富集分析

首先，加載DOSE包并使用該包中的基因集作為我們分析的對(duì)象：

library("DOSE")
data(geneList, package = "DOSE")
g_list <- names(geneList)[1:100]
head(g_list)

利用enrichKEGG函數(shù)進(jìn)行富集分析，對(duì)指定的基因列表進(jìn)行KEGG通路富集分析，并設(shè)置p-value和q-value的閾值。

library("clusterProfiler")
library("org.Hs.eg.db")

ek <- enrichKEGG(gene = g_list,
?? ? ? ? ? ? ? ? organism = "hsa",
?? ? ? ? ? ? ? ? pvalueCutoff = 0.05,
?? ? ? ? ? ? ? ? qvalueCutoff = 0.05)
write.table(ek, file = "ek.txt", sep = "\t", quote = FALSE, row.names = FALSE)

?

注意：

·?pvalueCutoff = 0.05：設(shè)置p-value的閾值，用于篩選富集分析結(jié)果中統(tǒng)計(jì)顯著的通路。在這里，設(shè)定為0.05，即p-value小于0.05的通路被認(rèn)為是顯著富集的通路。

·?qvalueCutoff = 0.05：設(shè)置q-value的閾值，q-value是經(jīng)過(guò)FDR校正的p-value，用于控制多重假設(shè)檢驗(yàn)的錯(cuò)誤發(fā)現(xiàn)率。同樣，設(shè)定為0.05，即q-value小于0.05的通路被認(rèn)為是顯著富集的通路。

至此我們已經(jīng)完成了富集分析的流程，讓我們一起來(lái)看一下分析的結(jié)果吧！

需要注意的是結(jié)果中g(shù)eneID為EntrezID。?

head(ek@result$geneID)

將EntrezID轉(zhuǎn)換為Symbol

使用setReadable函數(shù)將EntrezID轉(zhuǎn)換為Symbol：

ek2 <- setReadable(ek,
?? ? ? ? ? ? ? ? ? OrgDb = "org.Hs.eg.db",
?? ? ? ? ? ? ? ? ? keyType = "ENTREZID")

ggplot2可視化結(jié)果

首先，我們需要通過(guò)安裝tidyverse包并讀取保存的富集分析結(jié)果文本文件，對(duì)表格進(jìn)行處理和計(jì)算Enrichment Factor，在這里小果用tidyverse包的separate函數(shù)將GeneRatio和BgRatio的分子分母先分開(kāi)，再進(jìn)行計(jì)算：

library(tidyverse)
ek.rt <- read.table("ek.txt", header = TRUE, sep = "\t", quote = "") ?
ek.rt <- separate(data = ek.rt, col = GeneRatio, into = c("GR1", "GR2"), sep = "/")
ek.rt <- separate(data = ek.rt, col = BgRatio, into = c("BR1", "BR2"), sep = "/")
ek.rt <- mutate(ek.rt, enrichment_factor = (as.numeric(GR1) / as.numeric(GR2)) / (as.numeric(BR1) / as.numeric(BR2)))

篩選出前10個(gè)通路進(jìn)行可視化：

ek.rt10 <- ek.rt %>% filter(row_number() >= 1, row_number() <= 10)

ggplot2繪制dotplot

使用ggplot2包繪制Enrichment Factor與KEGG term之間的關(guān)系圖：

library("ggplot2")

p <- ggplot(ek.rt10, aes(enrichment_factor, fct_reorder(factor(Description), enrichment_factor))) +
??geom_point(aes(size = Count, color = -1 * log10(pvalue))) +
??scale_color_gradient(low = "green", high = "red") +
??labs(color = expression(-log[10](p_value)), size = "Count",
?? ? ? x = "Enrichment Factor", y = "KEGG term", title = "KEGG enrichment") +
??theme_bw()
ggsave("er.rt10_KEGG.png", width = 6, height = 4)

我們一起來(lái)看看效果圖吧！

繪制KEGG通路網(wǎng)絡(luò)圖

使用enrichplot包的cnetplot函數(shù)繪制KEGG通路網(wǎng)絡(luò)圖，并將結(jié)果保存為PDF文件。

library("enrichplot")

pdf(file = "net-pathway.pdf", width = 8, height = 10)
cnetplot(ek2)
dev.off()

效果圖如下：

以上教程小云介紹了如何使用R語(yǔ)言中的DOSE和clusterProfiler包進(jìn)行基因集的KEGG通路富集分析，并通過(guò)ggplot2和enrichplot實(shí)現(xiàn)可視化化。通過(guò)這些代碼，你可以將自己的基因列表進(jìn)行富集分析，并進(jìn)一步將結(jié)果進(jìn)行可視化以幫助解釋和理解數(shù)據(jù)哦！怎么樣，你學(xué)會(huì)了嘛？