微生物群落基於KEGG預測功能的豐度分布圖繪製

功能預測

在基於16S rRNA擴增子的研究中，因為了解微生物群落功能的需要，我們經常會使用一些軟體工具對16S rRNA擴增子數據進行功能預測，比如PICRUSt2。

這些功能預測的結果大多是基於KEGG資料庫，KEGG是一個整合了基因組、化學和系統功能信息的資料庫。

把從已經完整測序的基因組中得到的基因與更高級別的細胞、物種和生態系統水平的系統功能關聯起來是KEGG資料庫的特色之一。

與其他資料庫相比，KEGG 的一個顯著特點就是具有強大的圖形功能，它利用圖形而不是繁縟的文字來介紹眾多的代謝途徑以及各途徑之間的關係，這樣可以使研究者能夠對其關注的代謝途徑有更直觀全面的了解。

KEGG GENES資料庫提供關於在基因組計劃中發現的基因和蛋白質的序列信息。

KEGG PATHWAY資料庫包括各種代謝通路、合成通路、膜轉運、信號傳遞、細胞周期以及疾病相關通路等，此外還收集了各種化學分子、酶及以及酶促反應等相關信息。

KEGG Module資料庫是 KEGG 收集的一系列的功能單元，用於基因組注釋和生物學解釋。

KEGG Orthology (KO)系統通過把分子網絡的相關信息連接到基因組中，提供了跨物種注釋流程。

ko：表示通路，這個通路是不分物種的，相當於所有物種某一功能步驟的併集。

KEGG 資料庫包含 8 類生物代謝通路的分析：

• Global Map

• Metabolism, Genetic Information Processing

• Environmental Information Processing

• Cellular Processes

• Organismal Systems

• Human Diseases

• Drug Development

在每一個大類的代謝通路中，被系統分類為二、三、四層。

第二層目前包括有 39 種子功能；第三層即為其代謝通路圖；第四層為每個代謝通路圖的具體注釋信息。

豐度分布圖繪製

使用PICRUSt2等功能預測工具得到的KEGG路徑預測結果同樣是有層級的，當我們想要知道樣本中微生物群落的功能主要涉及哪些路徑時，就需要對這些具有層級的功能預測結果進行繪圖。

一般情況下在圖像中會展示第一和第二層級，因為三、四層級條目過多，在一幅圖中進行展示也看不清，通常會挑選感興趣的特定功能單獨呈現。

繪圖數據

繪圖數據是一個簡單的三列表格，包括第二層級的39個代謝路徑及其對應的第一層級信息，還有就是每個代謝路徑在所有樣本中的平均相對豐度。

⚠️因為目的是獲得所有樣本中的主要功能，因此繪圖數據需要提前計算一下每一個代謝路徑在所有樣本中的平均豐度。

圖像繪製

使用ggplot2的分面模式來展示不同代謝路徑所屬的第一層KEGG信息，應用條形圖表示每一個代謝路徑的相對豐度。

首先導入分析數據。

KEGG <- read.table("L2.txt",header = TRUE,sep = "\t")

這裡有一點需要說明，因為KEGG第一層分類中有幾個分類的名字特別長，因此需要處理一下，使其能夠實現自動換行。

library(ggplot2)library(reshape2)
swr = function(string, nwrap = 12){
paste(strwrap(string,width = nwrap),collapse = "\n")
}
swr = Vectorize(swr)
KEGG$L1 <- swr(KEGG$L1)

接下來就可以畫圖了。

p <- ggplot(KEGG,aes(Abundance,L2)) +
geom_bar(aes(fill = L1),stat = "identity",width = 0.6) +
xlab("Relative abundance (%)") +
ylab("KEGG Pathway") +
theme(panel.background = element_rect(fill = "white",colour='black'),
panel.grid.major = element_line(color = "grey",linetype = "dotted",size = 0.3),
panel.grid.minor = element_line(color = "grey",linetype = "dotted",size = 0.3),
axis.ticks.length = unit(0.4,"lines"),
axis.ticks = element_line(color='black'),
axis.line = element_line(colour = "black"),
axis.title.x=element_text(colour='black', size=8,face = "bold"),
axis.title.y=element_text(colour='black', size=8),
axis.text.x=element_text(colour='black',size=8,
axis.text.y = element_text(color = "black",size =8),
legend.position = "none",
strip.text.y = element_text(angle =0,size =12,face = "bold")) +
facet_grid(L1~.,space = "free_y",scales = "free_y")

保存圖片就OK了。

ggsave("L2.pdf", p, width = 183, height = 240, units = "mm")

png(filename="L2.png",res=600,height=7200,width=6000)
p
dev.off


PS：有些朋友可能想要把分面的文字放在圖像的左側，使其與代謝路徑的名稱放在一起，但是這個直接用ggplot2可能不行，因此把分面名稱放在左側之後，對其進行文字調整的參數就失效了，之前在群裡也討論過這個問題，據說是ggplot2的一個小bug，以後的更新中可能會修復。

後臺回復「功能預測」獲取示例文件和完整代碼。

10000+：菌群分析 寶寶與貓狗 梅毒狂想曲 提DNA發Nature Cell專刊 腸道指揮大腦

系列教程：微生物組入門 Biostar 微生物組 宏基因組

專業技能：學術圖表 高分文章 生信寶典 不可或缺的人