最近做細(xì)菌比較基因組學(xué)分析嗽桩，跑了一下Roary軟件，出來(lái)了一些結(jié)果和圖，看了看Roary軟件github上的作圖代碼机杜，使用的python語(yǔ)言對(duì)系統(tǒng)進(jìn)化樹(shù)和gene presence/absence數(shù)據(jù)的可視化，想要調(diào)整圖片的一些細(xì)節(jié)逼著自己去看python代碼衅谷，谷歌了 Biopython中的phylo模塊的使用文檔椒拗，依然是懵逼的，我想讓系統(tǒng)進(jìn)化樹(shù)顯示support值获黔，顯示并調(diào)節(jié)樹(shù)上的菌株名字的大小蚀苛，還是決定用自己稍微熟悉的R語(yǔ)言去解決問(wèn)題。
就是類(lèi)似這張圖要修改玷氏。
發(fā)現(xiàn)Y叔的ggtree可以實(shí)現(xiàn)類(lèi)似的圖堵未。
但我學(xué)習(xí)的時(shí)候沒(méi)有找到這個(gè)圖的完整數(shù)據(jù)，搞不清數(shù)據(jù)結(jié)構(gòu)是怎么樣的盏触，所以我決定用ggtree展示系統(tǒng)進(jìn)化樹(shù)渗蟹，ggplot2展示矩陣數(shù)據(jù)，然后使用patchwork兩個(gè)圖拼在一起赞辩。
1. 系統(tǒng)進(jìn)化樹(shù)可視化
#########系統(tǒng)進(jìn)化樹(shù)可視化#######
rm(list = ls())
getwd()
setwd("C:\\Users\\Administrator\\Desktop\\R_work\\07_BD177_genome_work")
list.files()
#載入包
library("ggtree")
library("treeio")
library("tidyverse")
#樹(shù)文件導(dǎo)入
BD177_tree <- read.newick("accessory_binary_genes.fa.newick",node.label = "support")
#系統(tǒng)進(jìn)化樹(shù)可視化
p1 <- ggtree(BD177_tree, branch.length = "none") +
geom_tiplab(size=2, align=TRUE, linesize=.5) +
xlim(0,30) +
geom_nodelab(aes(subset = !isTip,label=round((node/232), 2)),hjust=0.5,color="black", size = 2.5)
p1
2. gene presence/absence矩陣數(shù)據(jù)的可視化
#矩陣數(shù)據(jù)導(dǎo)入
gene_pre_ab <- read.csv("gene_presence_absence_heatmap_tidydata.csv", sep = ",")
#寬數(shù)據(jù)變長(zhǎng)數(shù)據(jù)
gene_pre_ab_long <- gene_pre_ab %>%
pivot_longer(cols = -strain, names_to = "gene", values_to = "value")
#將ggtree生成的系統(tǒng)進(jìn)化樹(shù)所有信息存入
p1b = ggplot_build(p1)
head(p1b$data)
#將gene_pre_ab_long數(shù)據(jù)中的strain與系統(tǒng)進(jìn)化樹(shù)數(shù)據(jù)中的label進(jìn)行一一對(duì)應(yīng)
gene_pre_ab_long = gene_pre_ab_long %>%
mutate(strain = factor(strain, levels=p1b$data[[3]] %>% arrange(y) %>% pull(label)))
#矩陣可視化
p2 <- ggplot(gene_pre_ab_long, aes(x = fct_reorder(gene, value,.desc = T), y = strain)) +
geom_raster(aes(fill=value)) +
scale_fill_gradient(low="#f6fafe", high = "dodgerblue3") +
labs(x="", y = "", title="") +
theme_bw() +
theme(axis.text.x=element_blank(),
axis.text.y=element_blank(),
plot.title=element_text(size=5)) +
guides(fill = FALSE )
p2
3. 拼圖
#拼圖與導(dǎo)出
library("patchwork")
library("export")
p1 + p2 + plot_layout(widths = c(1.5, 1))
graph2ppt(file="BD177_tree_heatmap.pptx", width=7, height=10)
Y叔的ggtree包還需要進(jìn)一步探索雌芽，看ggtree的幫助文檔發(fā)現(xiàn)功能十分強(qiáng)大。
